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ABSTRACT 

Observations of the intergalactic medium (IGM) suggest that quasars reionize Hell in the IGM at 
z«3. We have run a set of 190 and 430 comoving Mpc simulations of Hell being reionized by quasars 
to develop an understanding of the nature of Hell reionization and its potential impact on observables. 
We find that Hell reionization heats regions in the IGM by as much as 25, 000 K above the temperature 
that is expected otherwise, with the volume-averaged temperature increasing by ~ 12, 000 K and 
with large temperature fluctuations on ~ 50 Mpc scales. Much of this heating occurs far from 
quasars by photons with long mean free paths. We find a temperature-density equation of state of 
7— 1 « 0.3 during Hell reionization, but with a wide dispersion in this relation having ut ~ 10 4 K. Hell 
reionization by the observed population of quasars cannot produce an inverted relation (7 — 1 < 0). 
Our simulations are consistent with the observed evolution in the mean transmission of the Hell Lya 
forest. We argue that the heat input from Hell reionization is unable to cause the observed depression 
at z w 3.2 in the HI Lya forest opacity as has been suggested. We investigate how uncertainties in 
the properties of QSOs and of Hell Lyman-limit systems influence our predictions. 
Subject headings: cosmology: theory - intergalactic medium 



1. INTRODUCTION 

In the standard picture for the reionization history of 
the Universe, radiation from Population II stars ionized 
the intergalactic HI at z > 6 as well as the Hel, convert- 
ing the vast majority of the intergalactic helium to Hell. 
However, these stars cannot ionize Hell, and at z ~ 3 
quasars, with their harder UV spectrum, doubly ionize 
the intergalactic helium. To test this model, many ob- 
servations are tar geting high redshifts to probe hydrogen 
reionization (e.g.lFan et al.l 2002; T aniguchi et al.ll2005l: 
Bouwens et al.ll20q7t: iKashikawa et alJl2006t IStark et all 
2007t iTotani et al.l 120061 ). In this picture, Hell reion- 
ization occurs at redshifts for which there is substan- 
tially more data on the state of the intergalactic medium 
(IGM). 

In fact, a number of observations suggest that Hell 
reionization happened at z ~ 3. Two measurements 
of the mean transmission in the HI L ya forest have 
noted an upward bump at z ~ 3.2 (Bcrnar di et al 



120031: iFaucher-Giguere et al.ll2008aD . which lTheuns et al 
( 2002bf ) interpreted as arising from a temperature in 



crease of the IGM during H ell reionization (but see 
IFaucher-Giguere et al.l (|2008ah for alternative explana- 
tions). An increase in the average temperature of the 
IGM would also d ecrease the small-sc ale fl uctuations in 
the HI Lya forest. iRicotti et al.l (2000) and lSchaye et al.l 
(2000) measured the temperature from the widths of the 
narrowest lines in the HI Lya forest and claimed to have 
detected a sudden increase in the temperature of AT ~ 
10 4 K between z = 3.5 and 3. Photo-heating during Hell 
reionization is the only known process that could be re- 
sponsible for such an effect (e .g., lMiralda-Escude fc Reesl 
Il99l lAbeffc HaehneltlflQQff K 
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However, a subsequent study by 'McDonal d et al 
(1200 ll) using a similar method and Zaldar riaga et al. 
(2001) using the HI Lya forest power spectrum did 
not confirm this sudden increase in temperature, but 
rather found a constant temperature at mean density 
of T w 17,000 K for 2 < z < 4. Temperatures of 
17,000 K are difficult to explain without Hel l reioniza- 
tion occurring at z ~ 3 (|Hui fc Haiman| [2u03). and it is 
unclear whether as sudd en an increase in tem perature as 
IRicotti et ail ((2000) and lSchave et alJ (|2000f ) find is even 
expected theoretically. 

If a substantial fraction of the helium is in Hell (> 
1%), this would produce a Gunn-Peterson absorption 
trough in the spectra of high-redshift quasars at wave- 
lengths blueward of Hell Lya. Observations of Hell 
Lya forest absorption at 2.8 < z < 3.3 find 10s of 
como ving Mpc regions with no detected transmission 
(Jakobsen et alJll994HDavidsen et al.lll996HHogan et all 
119971 iReimers et al.1 Il997t iHeap et all |2000D . These 
troughs may signify the presence of diffuse intergalac- 
tic Hell. However, current data, which consist of only 
a few quasar sight- lines, do not rule out the intergalac- 
tic Hell being primarily ionized and in photo-ionization 
equil i brium with a weak background (jGiroux fc Shulll 
[J9971 [Fardal et alJll998t IHeap et al.ll200fl . The Cosmic 
Origins Spectrograph, which NASA plans to install on 
the Hubble Space Telescope in 2009, will increase the 
quantity and quality of Hell Lya forest sight-lines. 

As the intergalactic Hell becomes progressively more 
ionized, the extragalactic UV background will harden 
around the ionization energy of Hell at 54.4 eV. This 
hardening will affect the ionization state of intergalac- 
tic metals. ISongailal (J1998I ) observed a sharp evolu- 
tion at z rs 3 in the column density ratios in SIV 
(Ionizati on Potential = 45 . 1 eV) to CIV (64.5 eV) ab- 
sorbers. iBoksenberg et all (|2003|) found evidence for a 
more gradual hardening of the background between 2 < 
z < 4 from the column density ratios of NV (98 eV) to 



CIV. Finally, by simultaneously fitting to multiple metal 
lines that origin a te fro m th e same absorption syste ms, 
lAgafonova et al.l (|2005| ) and lAgafonova et all (|2007ft in- 
ferred a background spectrum that is hardening at z w 3 
near 4 Ry. 

Measurements from the z ~ 3 HI Lya forest ignore 
the effects of a patchy Hell reionization process. For 
example, estimates of the photo-ionizing background and 
the IGM temperature from the forest assume a power-law 
temperature-density relation. Cosmological parameter 
studies from the HI Lya forest power spectrum make 
a similar assumption. Different regions can have vastly 
different Hell reionization histories, resulting in a more 
complicated distribution of temperatures and pressure- 
smoothing scales than is commonly adopted. Realistic 
simulations of Hell reionization will help quantify the 
level at which this process biases these measurements. 3 

The aim of this paper is to run realistic simulations of 
Hell reionization to understand the morphology of this 
process as well as its effect on observables. We concen- 
trate primarily on its impact on quantities that are sensi- 
tive to the IGM temperature, but we also study the effect 
of Hell reionization on the transmission in the Hell Lya 
f orest. 

ISokasian et~aH (|2002f ) and iPaschos et all (|2007t ) have 
performed the most realistic simulations of patchy Hell 
reionization to date. There are several differences be- 
tween our work and these earlier investigations. Both of 
these studies employed volumes < 100 3 comoving Mpc 3 . 
Here, we examine Hell reionization in 186 3 and 429 3 
comoving Mpc 3 volumes, providing a more representa - 
tive cosmic sample . Howe ver, both lSokasian et all (2002) 
and IPaschos et all (|2007l ) simulated Hell reionization as 
a post-processing step on top of cosmological simulations 
that included gas dynamics. Our study instead uses N- 
body simulations, which result in a less realistic model for 
the gas distri bution, but afford a la rger dynamic range. 
Furthermore, ISokasian et al.l (|2002f ) assumed sharp ion- 
izing fr onts and ignor e d the detailed temperature evo- 
lution. IPaschos et al.l (|2007f) did not calculate the gas 
temperature self-consistently. Our calculations capture 
the width of the ionization fronts and the temperature 
in a consistent manner. Finally, in contrast to previous 
studies, our work presents a large set of radiative transfer 
simulations in order to survey the parameter space. 

In Section [21 we describe the details of our code. The 
models for the quasar sources are described in Section 
[3l Section [4] presents the simulations. Finally, Section [5] 
addresses the implications Hell reionization has on ob- 
servations of the HI and Hell Lya forests. 

Throughout, we use a ACDM cosmology with n s = 1, 
a s = 0.8, Q m = 0.27, fl A = 0.73, n b = 0.046, and 
h = 0.7, consistent wi th the most recent WMAP results 
(|Komatsu et alj I2008T) . All distances are in comoving 
units unless specified otherwise. An overbar over a vari- 
able signifies a volume average, and xy is the fraction of 
helium/hydrogen that is in ionization state Y . 

2. ALGORITHM 



Our work employs a new ray-tracing code that 
has been adapte d significantly from the one origi- 
nally presented in ISokasian et al.l (|2001l ) and refined in 
iMcQuinn et al.l (|2007f ). This code performs cosmological 
radiative transfer as a post-processing step on top of N- 
body or Smooth Particle Hydrodynamics density fields. 
In this study we use N-body simulations. The differences 
betwe en ou r code and thos e descr ibed in ISokasian et alj 
(|2001h and IMcQuinn et al.l (|2007h include that it tracks 
photons in multiple frequency bins, that it calculates the 
temperature of the gas, that it does not assume sharp 
ionizing fronts, and that it is parallelized over shared 
memory. This section describes the details of our code, 
and Appendix IBI discusses various tests of it. 

At the beginning of each timestep, the code grids the 
particles from an N-body snapshot, and inputs a list of 
source positions and luminosities. Then, the code adjusts 
the neutral fraction to account for the total number of 
recombinations that will occur over the ensuing timestep. 
Next, the code casts rays from every source, randomiz- 
ing the order of the sources in this loop as well as the 
direction in which the rays are cast. Every ray carries a 
set number of photons. Initially, N lay = 12 x 4 L rays are 
cast with L — 5 or 6 for an isotropically emitting source, 
with these rays uniformly pixelating the unit sphere using 
the HealPIX algorithm as described in lAbel fe; Wandera 
(2002). 4 Each source must wait for all of the other 
sources to send a ray before it sends a second ray. Rays 
are split adaptively as they traverse the box. At a min- 
imum, one ray from a source intersects every cell face 
for cells within the light travel time of a ray. (Although 
cells that are tens of Mpc from a quasar typically receive 
multiple rays from the source in a timestep, owing to the 
large numbers of rays that are initially sent.) 

A ray travels either until its photons are expended 5 
or until it has traveled a distance cAt, where At is the 
simulation timestep. Rays that have traveled c At are 
stored onto disk until the following timestep, at which 
time the stored rays are redshifted in frequency, ran- 
domly ordered, and recast. Rays that have traveled a 
total of 1.5 box lengths are terminated and their pho- 
tons are added to the background (Section 12.11 Most 
cosmological radiative transfer codes allow rays to travel 
further within a timestep than c At. To study Hell reion- 
ization, it is crucial to capture light-travel effects. 

For the calculations presented here, rays carry photons 
with energies between the ionization potential of Hell 
(E'Heii = 54.4 eV) and the energy that has a mean free 
path (m.f.p.) equal to 3.5 times the size of the 186 Mpc 
simulation box (1.8 times for the 429 Mpc box), which 
works out to E 1 ~ 500 eV in the 186 Mpc box at z — 3. 

4 While the value for iV ray may seem large, most of the compu- 
tation time is spent tracing rays that are far from the source so 
that the speed of the algorithm depends weakly on the number of 
rays that are initially cast. 

5 The termination criterion is that the number of photons in 
the ray be e times the number of atoms in a fully neutral cell 
with baryonic overdensity equal to zero, where we set e ~ 10 — 4 — 
10~ 5 with the exact value depending on the run in question. The 
photons from terminated rays are placed in the background. We 



have run a simulation with with e 



10" 



to test convergence 



3 Of note. iLaFc t al. (2006) used simple models for Hell reion- 
ization to show that it has a surprisingly small effect on the Lya 
forest power spectrum on large scales, modifying it at a < 5% level 
for wavevectors k < 5 comoving Mpc -1 . 



(simulation Lid in Table 1) and find no notable differences with 
simulation LI for which e = 8 X 10 — 4 . In fact, the temperature 
and ionization fields are almost visually indistinguishable between 
these two simulations. 



The number 3.5, while somewhat arbitrary, assures that 
the flux is uniform on the box scale for photon energies 
that are not tracked by the rays. Photons up to 5 keV 
are put into an ionizing background f£ )2.ip . Many of the 
background photons are never absorbed because photons 
with £ 7 w 0.8 x [(1 + z)/4] 1 /2 s ^ 3 n keV have a m.f.p. 
equal to the Hubble scale. 

Each ray carries some number of photons with en- 
ergy at Ei for 1 < i < n v . This study typically uses 
n v = 5. The value of E~ i of the photons within an en- 
ergy bin is specified to conserve photon number and to- 
tal energy given the source spectrum. The centers of 
the energy bins that are tracked by the rays are roughly 
(60, 70, 100, 140, 270) eV for our fiducial spectrum in the 
186 Mpc box. In Appendix [B| we show that 5 is the min- 
imum number of frequency bins that provides adequate 
convergence. The final volume-averaged temperature of 
a simulation with a much larger number of energy bins 
is within 300 K of the same run with n v = 5 energy bins, 
and the ionization and temperature fields are similar. 

Ideally, radiative transfer codes would employ a 
timestep that is much smaller than the photo-ionization 
equilibrium timcscalc. For cosmological ray tracing 
codes, such small timest eps can be c omputationally in- 
tractable (although, see iTrac fc Cenll2007f h Our algo- 
rithm is designed to operate with a timestep that can be 
much larger than the equilibration timescale; although, 
in this study it is typically used in the regime where 
the two are comparable. The timestep used in our sim- 
ulations is At w 10 Myr. Compare this to the photo- 
ionization equilibration timescale (rHeii + o^e) -1 , where 
ThcII is the Hell photo-ionization rate and a is the re- 
combination rate coefficient. This timescale, which in 
most cases is dominated by the THeii term, is observed 
to be 3 - 30 Myr at 2.5 < z < 3 (i.e. -15 < log 10 r Hc ii < 
— 14 in c.g.s.), but it can be much shorter near quasars. 
An accurate solution can be obtained with our time steps 
as demonstrated in Appendix [Cl 

Converging to the solution also requires capturing the 
correct order that rays intersect a cell. The number of 
absorptions a ray experiences depends on whether it ar- 
rives before or after previous rays that have intersected 
a cell. If the timestep is small enough, such that the ion- 
izing front does not move significantly over a timestep, 
the order rays hit a cell does not matter. Our code is 
not always operating in this limit. It attempts to reach 
a converged solution by having many rays from a source 
intersect a cell that is within ~ 30 Mpc from the QSO 
and by randomizing the order in which the rays are cast. 

2.1. Background 

Hell-ionizing photons with A m f p > 3.5 Lbox, where 
A m fp is calculated assuming a homogeneous IGM with 
£Hcii = 1, are put into a pool of background photons, 
which we specify with 10 frequency bins. 6 Five keV is 
the maximum photon energy that is included in the back- 
ground. Averaged over a timestep, each cell receives a 
background flux of 

/" = c 73 — — I 1 ~ ex P(- T ^)] . 00 

^box T " 

6 This number is for the 186 Mpc box. We use the criterion 
-Wp > 1.8 -Lbox for the 429 Mpc box. 



where A^k^ is the total number of background photons 
in frequency bin i/, Lbox is the size of the box, and t v is 
the average optical depth of a photon traveling a distance 
c At in the simulation volume. 7 



2.2. Temperature 

Our code tracks the temperature evolution of the 
gas, which is govern ed by the differential equation (e.g., 
iHui fc Gnedinlfl99l 



2T dS b T dJ^Xi 



— = -2HT 

dt 3 A 6 dt J2 X 



ill 



+ - 



dQ 



3k B n tot dt ' 
(3) 

where Xi is defined such that the number density in 
species i is (1 + $) Xi pb/m p , dQ /dt is the heating rate, 
and ritot is the total number of gas particles. 8 The first 
and second terms on the right-hand side of equation ([3]) 
account for the effect of adiabatic heating/cooling, the 
third term describes the change in the number of gas 
particles owing to ionizations and recombinations, and 
the final term accounts for radiative heating and cooling. 
The most important process that contributes to dQ/dt is 
photo-heating, but atomic cooling, recombination cool- 
ing, dielectric recombination cooling, collisional cooling, 

7 A separate background is also included to determine the 
ionization state of hydrogen. While we could include this ion- 
izing background self-consistently by extrapolating the 4 Ry 
luminosity of quasars in the box to 1 Ry, there is mounting 
evidence that there is a substantial stell ar contribution to 
the ionizing background at these redshift s |[Madau_et_alJ 119991; 
Stcidel ct al. 2001; Bokscnbcrg ct al. 2003; Sokasian ct al. 2003; 
Fauchcr-Gigucrc ct al. 2008b). Therefore, we instead adopt an em- 
pirical approach, utilizing observations of the HI photo-ionization 

rate (fni) fro m the H I Lya forest. 

Bolton ct al. (2005), [Becker et al.l <|2007I 1. and especially 
Fauchcr-Gigucrc et al.l (l2008bf) measured the hydrogen photo- 
ionization rate (rjji) to be flat as a function of redshift using the 
flux decrement method on 2 < z < 4 HI Lycr forest data and to 
within a factor of ~ 2 given by 



r m (z) = io- 



(2) 



We use this function for rjji(z) throughout this paper and for all 
simulated redshifts. The background level primarily affects our 
calculations regarding the HI Lyo forest opacity and does not in- 
fluence our conclusions. 

The approximation of a global background is excellent for track- 
ing the photo-ionization of HI at the redshifts relevant for 
Hell reionization because the m.f.p. for hydrogen ionizing pho- 
tons is measured to be l m f p = 300 [(1 + z)/4}~ 3 - 2 Mpc 
IjMeiksin &: W hite 2004), with present uncertainty at a factor of 
two level (Fauchcr-Gigucrc ct al. 2008b). A rough estimate is that 
there is 1 galaxy per Mpc 3 and one quasar per 30 3 Mpc 3 , such 
that there are thousands of sources that contribute to the back- 
ground within £ m f p . This implies that spatial fluctuations in the 
background are small at z < 5 (Croft 2004; McDonal d et al.D2005l) . 
For similar reasons, this approximation is also valid for the hard 
Hell-ionizing photons that are treated as a global background by 
our code. 

8 Several published studies that solve for T use an incorrect 
form for the third term on the R.H.S. of equation (|3j - a factor of 
—2/3 different from what appears in this equation. This mistake 
stems essentially from writing the first law of thermodynamics as 
dq p = de p + P d(l/n), where q p is the external heating per particle, 
e p is the energy density per particle, and P is the pressure. The 
correct form for this equation is dq p = de p + P d(l/(fj,m p n)), where 
fi is the mean molecular weight and m p is the proton mass. Unlike 
in the former equation, the second term in the correct version does 
not depend on the ionization state of the gas. This mistake results 
in the temperature of the gas increasing during ionization processes 
by too large an amount. 



and Compton cooling off of CMB photons are also in- 
cluded in our calculation. 9 We use an implicit solver to 
obtain a solution to equation ([3]). 

To evaluate the Lagrangian time derivatives in equa- 
tion J3J within a cell, the locations of particles at the 
beginning of the timestep are used to calculate the ini- 
tial values for the ionization, density, and temperature. 
Namely, each particle initially has the ionization, den- 
sity, and temperature of the cell that it was in at the end 
of the previous timestep. Summing up these previous 
values for all the particles within the cell (with the ap- 
propriate weighting to conserve energy) gives the initial 
values of these scalars. 

The density field is calculated from the gridded N- 
body particles. The dynamics of N-body particles are 
different from that of gas particles. Therefore, it is ques- 
tionable whether this algorithm will be able to predict 
the correct temperature. However. iHui fc Gnedinl (|1997f ) 
demonstrated that the temperature evolution produced 
from smoothing dark matter simulations at the Jeans 
scale provides good agreement with the temperature evo- 
lution seen in hydro dynamic simulations. In addition, 
IHui fc Gnedinl (|1997l ) showed that even using linear per- 
turbation theory to evolve the density field provides rea- 
sonable agreement with the evolution of the T-Aj, rela- 
tion in hydrodynamical simulations, suggesting that ex- 
actly capturing the Jeans scale is not crucial. A compar- 
ison | _of_Hie_tejn£erature_e volution our code predicts with 
the IHui fc Gnedinl (|1997f l analytic formula is presented 
in Appendix 151 

The simulations in our study capture the Jeans scale 
to varying degrees. [Although, more appropriate may be 
the filtering scale, w hich is typically a factor ~ 2 smaller 
(jGnedin fc Huilfl998h .] The Jeans mass is given by 



M 3 = 9.6 x 10 9 A 



-1/2 / T 
\ 10 4 K 



.3/2 



1 



-3/2 



Mr. 



(4) 

Compare this mass to the mass within a grid cell in our 
simulations: 

M cel] - 1.4 x 10 10 A 6 ( Lbox -P^ M Q . (5) 
V186Mpc N ccns J ° V ' 

These masses are comparable for Ab = 1 for our fiducial 

resolution and box size (£box = 186 Mpc and N cc \\ s = 

256). For the 512 3 runs in the 186 Mpc box, M ceU < 

Mj(z = 3) for A b = 1, and, in the 429 Mpc box, 512 3 

resolution achieves M ce ii « Mj(z — 3) for Aj = 1. Most 

important, our conclusions do not change if we use the 

512 3 mesh rather than 256 3 . 

2.3. Secondary Photo-ionizations and Heating 

We assume that the excess energy above -EhcII from 
a photon that photo-ionizes a Hell ion goes into heat- 
ing the gas. In reality, the secondary electron pro- 
duced by the photo-ionization may go on to colli- 
sionally ionize another atom, excite a bound elec- 
tron, or suffer coulomb co llisions that heat the gas 
(IShull fc van Steenberdfl985h . 

9 We assume a power-law of 1.6 in energy flux to calculate the HI 
photo-heating. We neglect Hel photo-heating in this calculation. 
The heat input from Hel photo-heating is » 25% that of the HI 
photo-heating. Uncertainties in the spectrum of the photo-ionizing 
radiation result in larger uncertainties in the photo-heating rate 
than the Hel contribution to it. 



However, IShull fc van Steenberd (|1985l ) find that 99% 
of the excess energy of a 3 keV X-ray photon that photo- 
ionizes a Hell ion goes into heating the gas for a situation 
in which xhi = 0.05, ShcI = 0.05, and xhcII = 0.95. This 
percentile will increase for the smaller Shi, a^Hei, and pho- 
ton energies that are most relevant during Hell reioniza- 
tion. Therefore, the assumption that all the excess en- 
ergy goes into heating t he ga s is excellent. Furthermore, 
IShull fc van Steenberd (1985) find that interactions with 
Hell are unimportant in their calculations. Most of the 
excess energy is converted into heat because collisions 
with electrons dominate the energy loss of the secondary 
particles once hydrogen is reionized. 

Finally, Compton cooling of the energetic electrons off 
of the CMB, which was not in cluded in the calculations of 
IShull fc van Steenberd (|1985l ). is a subdominant energy 
loss mechanism at z ~ 3 for the releva nt gas densities 
and photon energies (jMadau et al.ll2004} ). 

2.4. Heat Input Estimates 

Let us develop an understanding for how much heat- 
ing is expected from Hell reionization. An Hell ion- 
ization front heats u p the gas behind it by (e.g., 
I Abel fc HaehndtllT999l) 



LAT: 



2>Hc 



dE 



3(8-5!hc) V £hoII E 
(E-Eboi) cT He ii(.E) Ji(E) e- 

(f°° dE 
/ -ft OlMlCE) ME) e 



(6) 



where Ji is the unabsorbed intensity emitted by the 
source, Ji(Ey) exp(— te) is the incident spectrum, te is 
the optical depth from the source to the gas parcel, Yu e 
is the primordial helium mass fraction, and CHeii is the 
Hell photo-ionization cross section. 



For i£ C 1 and using that CHeii 
equation ([6]) becomes 



(E 1 — Eflell)' 



AT & - 



-EhcII 



3 27fc 6 (2 + a uv ) 



= 4500 



3.5 



«uv 



K, (7) 



where quv is the spectral index of Ji(E^), £hcIi/(2 + 
auv) is the average excess energy of a photon above 
Eu e u, and 27 is the number of particles over which this 
energy is distributed assuming Ihc = 0.25. If te is ap- 
preciable, equation (0) still has some relevance, except 
replace auv with the effective spectral index of the inci- 
dent spectrum, which will be harder, resulting in larger 
AT. 

The value 4500 K in equation ((7} is quite small, and it 
might be difficult to reconcile such a small AT with mea- 
surements of AT in the Lyman-a forest (assuming Hell 
reionization is the cause of the temperature increase) . It 
turns out that this number is an underestimate for the 
total heating during Hell reionization. If we assume all 
photons up to £ max are absorbed then the heat injection 
is 



AT w 31, 000 



/ 0.5 



\auv - 1 



1 



auv^ST 



E, 



auv-l 



K, (8) 



where we have kept only the leading order in -EhoIi/ E max . 
This equation can be derived by setting cthoII = 1 and 
te = in equation ([6]) and performing the integrals. The 



heating implied by equation {5} is much larger than by 
equation ((7j). Several hundred eV photons have optical 
depth unity on roughly the 186 Mpc box scale assuming a 
homogeneous IGM with ^hoII = 1, so these photons will 
be absorbed somewhere in the IGM if they are produced 
during Hell reionization (a photon travels 850 Mpc be- 
tween z = 4 and 3, roughly the interval over which Hell 
reionization occurs in our simulations). 

If we set -©max = 350 eV [which has r ~ 1 over 200 Mpc 
for ShcII = 1] yields AT = 13, 000 K using equation $E§ 
and auv — 1-5 (and 15, 000 K if we had not expanded in 
Eneii/E max ). This number is comparable to the average 
temperature increase seen in our simulations discussed 
in Section 2J 

2.5. Recombination Rates 

The recombination timescale for Helll is £ re c,Heiii = 
[«a(T) n e (z)}^ 1 , where «a is the relevant Case A re- 
combination coefficient and n e is the mean electron den- 
sity. If t r0 c,Hciii is expressed in terms of the age of the 
Universe, t un i « 2/3 H(z)^ 1 , it yields 



ircc.HcIII 



0.6 



T 



10 4 K 



0.7 



1 



-3/2 



A," 1 , 



(9) 



where A;, is one plus 5b - the overdensity in gas - and 
we have assumed that the intergalactic hydrogen is fully 
ionized. The Hell recombination timescale is 5.4 times 
shorter than this timescale for hydrogen at T — 2 x 10 4 K. 
This fact is important because additional recombinations 
require more ionizations to reionize the IGM, which re- 
quires more energy injection. Although, the recombina- 
tion timescale is comparable to the gas cooling timescale, 
so that the additional heating does not lead to signifi- 
cantly larger t emperatures in regio ns that recombine and 
are reionized (jBolton et al.ll2008bh . 

Since the number of HcIII recombinations during Hell 
reionization is large, recombination radiation could con- 
tribute significantly to ThcII- 10 Unfortunately, HcII- 
ionizing photons produced by recombinations to the 
ground state of Hell are no t followed by our radiative 
transfer code. iFardal et al.1 (j!998ft showed that this ra- 
diation could contribute to the photo-ionization rate at 
the ~ 20% level. Tracking this radiation with our code 
would be prohibitively expensive because it would re- 
quire treating all cells as sources of ionizing radiation. 11 
Instead, our code either assumes that this recombina- 
tion radiation contributes locally to ionizations, which it 
does by using the Case B recombination coefficient for 
Hell, or that these photons do not contribute to ion- 
izations of the diffuse IGM, which it does by using the 
Case A recombination coefficient. The code also uses the 
corresponding recombination cooling rates for these two 

10 Ground state recombination radiation results in a line pro- 
file of the form exp[— (E 1 — EY{ c ii)/kT] for E 1 > -EhcII, such 
that most photons have excess energies of ~ 1 eV above -EhoII 
for characteristic temperatures of T ~ 10 4 K. A photon will travel 
80 [(1 + z)/^" 1 / 2 ([Ej - -E H oll]/leV) Mpc prior to redshifting 
below -EhoII- This distance is longer than the m.f.p. for such pho- 
tons to be absorbed, and, therefore, they will typically be absorbed 
either in dense clouds or by diffuse intergalactic Hell. 

11 T wo recent ray tracing codes outlined in Pawlik & Schavel 
(|2008I 1 and I AltaT ct al. (2008) have developed techniques to allevi- 
ate this problem. 



cases. Treating Hell recombinations with the Case A co- 
efficient results in a hotter IGM; more ionizing photons 
from QSOs are required to ionize the IGM, and the QSO 
photons are harder than the photons from recombina- 
tions to the ground state. Neither Case A nor Case B 
is correct in detail, and the most appropriate choice for 
the recombination coefficient depends on the ionization 
state of the gas. 

During Hell reionization, a large fraction of the ion- 
izing photons from recombinations will ionize the dif- 
fuse IGM since the Hell bubble size is comparable to the 
m.f.p. to be absorbed in dense systems. In this case, 
Case B is the better choice. Since this paper is most in- 
terested in studying the Hell reionization process, most 
of our simulations use the Case B recombination coeffi- 
cient. However, we have run a case that uses the Case A 
coefficient (simulation D4 in Table 1). This simulation 
results in a slightly higher volume-averaged temperature 
at the end of Hell reionization, with the difference being 
AT ps 1000 K when compared to a similar simulation 
that uses the Case B rate. The character of the ioniza- 
tion and temperature fields are similar between these two 
simulations. 

2.6. Mean Free Path of Ionizing Photons 

At our grid scale, subgrid fluctuations may play an 
important role in absorbing and hardening the typical 
spectrum - "filtering" - of the Hell-ionizing radiation. In 
Appendix [Al we investigate how dense clumps filter the 
radiation, and we quantify how well this effect is captured 
in our simulations. Systems that have HI column den- 
sities Ahi ~ 10 15 cm -2 are responsible for filtering Hell 
ionizing radiation in Helll regions during Hell reioniza- 
tion. These systems are much less overdense and more 
diffuse than those that limit the m.f.p. of HI Lyman- 
limit photons , having overdensities of 5b ~ 10 at z ~ 3 
(Schayc 2001), and, therefore, they can be captured with 
lower resolution simulations than their counterparts for 
HI Lyman-limit photons, which have 5b ~ 100. In Ap- 
pendix [D] we show that our gridded density field has 
some success reproducing the column density distribu- 
tion of these systems when compared to high-resolution 
hydrodynamic simulations. However, even though the 
column density distribution is reproduced, the bias and 
density of these systems is altered compared to high res- 
olution hydrodynamic simulations. In addition, our ra- 
diative transfer code may systematically over-ionize these 
systems, as described in Appendix[D] Therefore, we sup- 
plement our calculations with two prescriptions (Meth- 
ods A and B) to study the effect that dense systems have 
on filtering the ionizing radiation. These methods are de- 
scribed in detail in Appendix lAl and briefly here. 

Method A assumes that the high-column density Hell 
systems that are resolved in our simulations are in 
photo-ionization equilibrium. We argue that this ap- 
p roach is reasonab l e in A ppendix lAl Method B uses 
a lHaardt fc Madaul (|1996l )-like model that takes the dis- 
tribution of iVm measured from Lyman-a forest spectra 
and, given a model for the density of these systems and 
the local value of the photo-ionization rates in the sim- 
ulation, infers the Hell column density distribution and, 
therefore, the opacity in these systems. Further, this 
model assumes that these absorbers are associated with 
halos in our simulation that have m > 7 x 10 9 Af n and 



() 



places the absorbers at the locations of these halos in 
the simulation box. The reader might be wary of such 
methods for supplementing the resolution, but we find 
that these prescriptions have only a minor impact on the 
final results (Section 14. 3p . 

If much of the intergalactic helium is in Hell, the m.f.p. 
is often limited by diffuse gas rather than by dense sys- 
tems. In the limit that diffuse regions near the mean 
density dominate the opacity, the m.f.p for a photon with 
energy E 7 is given by 



A 



mfp 



5 X 



I *^1 



HeI1 VlOOeV 



1 + z 



Mpc. (10) 



The value of A m f p scales strongly with E 7 in this limit. 
A photon with E-y — 55 eV has A m f P = 0.8 Mpc at z = 3, 
whereas one with E~ f — 200 eV has A m f p = 40 Mpc. The 
average energy of a photon for our fiducial UV power- 
law index of auv = 1-6 is 150 eV. However, half of the 
Hell-ionizing photons have 54.4 < E 1 < 84 eV. 

3. IONIZING SOURCE MODEL 

Quasars are the leading candidates for reionizing Hell. 
Cooling radiatio n from massive ha los is the next most 
likely contender (|Miniati et al.ll2004r ). However, the large 
fluctuations in the Hell to HI column density ratios at 
z ~ 3 strongly disfavor cooling radiation as being the 
dominant source of Hell-ionizing photons. Too many 
sources of cooling radiation are present wit hin one mean 
free p ath to source these large fluctuations (|Bolton et al.l 
2006). Finally, there are other more exotic sources that 
could potentially ionize Hell at z < 6. Possibilities in- 
clude an unknown source of > 4 Ry radiation that is 
generating the Hell 1640 A recombination line observed 
in t he composite spectrum of z ~ 3 Lyman bre ak galax- 
ies (|Shaplev et alj|2003t iFurlanetto fc Ohll2008t) 12 . or the 
bi-products of decaying/annihilating light dark matter 
particles. These more exotic scenarios are probably also 
in conflict with observations for the same reason that 
cooling radiation is disfavored. Our study concentrates 
on Hell reionization by quasars. 

Detailed observations of quasars in the past couple 
decades have provided many clues into the nature of 
these extremely luminous objects. We now know that 
quasars are powered by accretion onto super-massive 
black holes. Presently, the z < 3 luminosity function 
of QSOs is well constrained over a few decades in lumi- 
nosity in both the optical and the X-ray. The clustering 
of these objects reveals that optically selected quasars 
reside in h alos with m ~ 10 12 M&, in depend ent of their 
luminosity (|Porciani et al.ll200J iCroom et al.| :2005). Al- 
though, clustering is currently measured over only about 
a decade in luminosity. For the analysis in this study, 
the biggest uncertainties in modeling QSOs are their life- 
times, tqsOj the fraction of Hell ionizing photons that 
escape the environs of a QSO into the IGM as a function 
of direction, and their spectrum between 4 Ry and 1 keV. 

Figure [T] plots the bolometric luminosity function of 
quasars for four redshifts. The black, solid curves are the 



12 IFurlanetto fc Ohl ([20081 find that only if Hell ionizing photons 
can escape with / esc Hell ^ 0-5 from galaxies can they contribute 
a significant fraction of the ionizing background. These values 
are much larg er than the theoreti cal expectation for / osc jjoll from 
galaxies (e.g.. [Gnedin et al.H2008l) , 
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Fig. 1. — Comparison of models and measurements of the bolo- 
metric QSO luminosity function at 4 redshifts. The data points 
are from a compilati on of observations across different wavebands 
(Hopkins et al. 2007), using a spectral template to infer the bolo- 
metric luminosity function. The solid curves are the lHopkins et ahl 
1120071) best-fit luminosity function (which the quasars in Method I 
are drawn from) and the dashed curves are the luminosity functions 
generated with Method II. The solid horizontal lines represent the 
number densities at which 1 QSO is in the 186 and 429 Mpc boxes. 



iHopkins et all (|2007l) best-fit model. 13 The cyan points 
with error-bars are infrared observations, the green are 
optical, the blue are soft X-ray, and the red are hard X- 
ray. The references for these observations are given in 
IHopkins et al.l (|2007f ). and the symbols follow the same 
conventions as described in their Table 1. Observations 
in different bands are included in this figure by assum- 
ing an intrinsic spectral template and an average column 
density distribution for the obscuring matter that is mo- 
tivated by X-ray measurements of Ahi- At z > 3, the 
luminosity function falls in amplitude. At z > 5 there 
is not much data to constrain the models. The abun- 
dance of QSOs between 3 < z < 4.5 shapes the duration 
and morphology of reionization in our simulations. In- 
sufficient HcII-ionizing photons are produced at higher 
redshifts to substantially ionize the Hell with the best- 
fit luminosity function model. In this model, ~ 2 Hell 
ionizing photons per helium atom are produced by z = 3, 
about enough to reionize the Hell. 

The faint-end slope of the the luminosity function 
dcf>/d\ogL is constrained to be rather flat at z « 3, with 
power-law index —0.5. This flatness means that most 
of the Hell-ionizing photons are produced by L ~ L* 
quasars. The horizontal lines in the top panel in Figure 
Q]demark the number densities that yield 1 quasar in the 
186 and 429 Mpc simulation boxes. Our simulations are 
large enough to contain many L* QSOs. 

Quasars are placed in our simulations with one of two 
methods: 

3.1. QSO Method I 
The lifetime of all quasars is fixed to be tqso, and 
quasars are placed in halos of mass m oc L 4 ' 3 . This 



13 The Hopkins et al. (2007) luminosity function is derived by 
simultaneously fitting to the infrared, optical, and X-ray luminosity 
functions. These fits model the observed quasar spectra with a 
bolometric luminosity-dependent intrinsic spectral template plus a 
model for the column-density distributions in TVjji and dust. This 
procedure re sults in a good fit to the observed luminosity function 
in all bands (IHopkins et al.H2007r ). 



scaling assumes QSOs shine at their Eddington Lumi- 
nosity (X c dd oc Mbh) with some duty cycle and that 
M b h ~ ct 4 ~ m 4 / 3 . 14 The number of QSOs in the box at 
a given luminosity is drawn from a Poisson distribution 
with m ean set by the o bserv ed luminosity function. We 
use the iHopkins et al.l ([2007) luminosity function evalu- 
ated at 912 A and extrapolate across the relevant band 
with a power-law auv- 

3.2. QSO Method II 

Recently, simulations of merging galaxies have been 
used to predict quasar light curves (Hop kins et alj 2005, 
2006), and these predictions have been remarkably suc- 
cessful at reconciling observations of the quasar luminos- 
ity function, quas ar clustering, and the unresolved X-ray 
background (e.g., Hopk ins et al.l [20061 ) . The bolometric 
light curve of the quasars in merger simulations can be 
roughly parametrized as 



Lboi(t) = 0.8 £ c dd(Vttbh) x 



rexp[i/r s ] t<0 

1(1 + t/r s )- b t>0 



(11) 



where mbh is the super-massive black hole mass, b = 
1.7 + 0.7 x log w (L cdd /10 12 L Q ), L cdd is the Eddington 
luminosity, and ts = 40 Myr. To map from Lboi(i) to the 
intrinsic luminosity a t the HI Lyman- limit, we take the 
I Hopkins et al.l (2007) iboi-dependent spectral template. 
Next, equation (filj) plus a relationship between halo 
mass and L edd (mbh) allows one to deconvolve from the 
observed QSO luminosity function the triggering rate per 
unit volume of quasars as a function of time. The trig- 
gering rate is defined as the rate quasars shine at their 
peak luminosity. We use the Magorrian relationship to 
relate a ga laxy's stellar mass to its super-massive black 
hole mass (Hopkin s et al.ll2007allbT ), and the halo o ccupa - 
tion model of iTinker et all (|2005h : IConroy_et_aTl (l2006h: 
IVale fc Ostrikerl (|2006l ) (for details, see Hopkin s~et al.l 
(|2008af) ) to map from halo mass to stellar mass. 15 Figure 
[jshow the luminosity function that this model generates 
(red, dashed curves). While it differs somewhat from the 
best-fit luminosity function (black, solid curves), it agrees 
well with the observations. To implement this model in 
our simulations, we switch on quasars in the simulation 
box using the triggering rate per unit halo mass described 
above, accounting for Poisson scatter. 

3.3. Discussion of QSO Models 

Figure [5] shows the Helll bubble size probability distri- 
bution function (PDF) in Method I and Method II. These 
bubble sizes correspond to the size of a bubble if all the 
QSOs sit in separate spherical Helll regions and there 
are no recombinations. All curves are calculated assum- 
ing auv = 1-6. For the thin curves, the PDF is volume- 
weighted, and, for the thick curves in the bottom panel, 



The specific mapping 



7.9 X 



10 12 [(5vL„)/10 46 ergs- 1 ] 3 / 4 [(1 + z)/4] x - s M , where vL v 
is evaluated at 1 Ry. Our results depend weakly on the mapping 
between L and m because the morphology of Hell reionization 
turns out to be determined primarily by Poisson fluctuations in 
the abundance of quasars rather than by their clustering. 

15 As a final ingredient, since there is scatter in the relation- 
ship between m^h an d Af stG n aI ., we assume that the actual mbh is 
log-normal function with dispersion 0.3 dex and logarithmic mean 
given by the Magorrian relationship. This step is necessary in order 
to reproduce the bright end of the luminosity function. 
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Fig. 2. — Helll bubble size distributions in our two QSO models 
assuming ojtjv = 1-6, isotropic emission, 1 QSO per bubble, and 
no recombinations. The black solid curves correspond to z = 2, the 
green dashed to z = 4, and the blue dotted to z = 6. Top panel: 
Thin curves are the volume-weighted probability distribution of 
bubble radii in Method I. These curves assume tqso = 40 Myr. 
Bottom panel: Thick curves are the volume-weighted probability 
distribution in Method II, and the thin curves are the number- 
weighted probability distribution. 

it is number-weighted. Method I has an infinite num- 
ber of quasars with L < L* and so the number-weighted 
PDF is not well-defined. 

The Helll bubbles from the quasars in Method II have 
a characteristic size of 35 Mpc in the volume-weighted 
PDF, and the probability distribution of Rb does not vary 
significantly with redshift. With Method I, the charac- 
teristic radius is a bit smaller, 20 Mpc for tqso = 40 Myr, 
and the dispersion in Rb is larger. Method I has a non- 
negligible probability for extremely large bubble sizes to 
exist, bubbles with Rb > 70 Mpc. These large bubbles 
are averted in Method II because brighter quasars have 
the shorter lifetimes in this method. The properties of 
Hell reionization do not change substantially if we em- 
ploy Method I or Method II in our simulations. 

Given the large observational uncertainty in the mean 
and distribution of auv (see Appendix |Ej), we leave these 
quantities as free parameters in both QSO Models. For 
each quasar in the simulation box, auv is drawn from 
a Gaussian distribution with mean auv and s.d. a a . 
Our fiducial model assumes auv = 1.6 and a a = 0.2 , 
whi ch is motivated by t he results of iTelfer et al.l (|2002l ) 
and lSteffen et all (|2006[ ) [Appendix [E]. We find that the 
amount of heating from Hell reionization is sensitive to 
the spectrum of QSOs. 

We also include a parameter £, a suppression factor for 
the number of Hell-ionizing photons that are emitted by 
QSOs. This parameter gives us a knob to tune in or- 
der to have Hell reionization occur at a desired redshift, 
and, given current uncertainties in the QSO luminosity 
function, can be adjusted at a factor of 2 level. 

4. SIMULATIONS 

We have run a series of simulations to investigate the 
morphology of Hell reionization. All of our simulations 
track quasars to Lboi > 10 43 or 10 44 erg s _1 , approxi- 
mately 3 decades below L*. Table Q] lists the simulations 



TABLE 1 
HeII Reionization simulations. 



Sim. a 


ibox (Mpc) 


LLS 6 


QSO Model 


Emission c 


c 


Grid 


Comments 


Dl 


186 




II 


iso 


0.75 


256 




D2 


186 




II 


iso 


0.75 


256 


no density fluctuations (Si, = 0) 


D3 


186 




11 


iso 


0.75 


256 


twice as many frequency bins as Dl 


D4 


186 




II 


iso 


0.75 


256 


case A recombinations for Hell 


LI 


186 


A 


II 


iso 


0.75 


256 




Lib 


186 


A 


II 


iso 


0.75 


256 


7 = at z = 6 


Lie 


186 


A 


II 


iso 


1.5 


256 


slightly earlier reionization than LI 


Lid 


186 


A 


11 


iso 


0.75 


256 


extremely conservative photon termination 


L2 


186 


A 


11 


iso 


1 


512 


high resolution run 


L3 


186 


B 


II 


iso 


0.75 


256 


includes A H i > 10 14 - 5 cm" 2 


XI 


186 


A 


II 


iso 


0.7 


256 


omits photons in background 


SI 


186 


A 


II 


beamed 


0.75 


256 


beam uses Gilli et al. (2007) model 


S2 


186 


A 


I 


iso 


2 


256 


TQSO = 10 Myr 


S3 


186 


A 


11 


iso 


0.75 


256 


assumes qtjv = 1.2 


Si 


186 


A 


11 


iso 


0.75 


256 


assumes ayv = 0.6 


S4b 


186 


A 


II 


iso 


1.25 


256 


assumes otjv = 0-6, no background 


Bl 


429 




II 


iso 


1 


256 




B2 


429 


A 


II 


iso 


1 


512 





a Unless specified otherwise, the initial conditions for these simulations are 7 — 1 = 0.3 and To = 10 4 K at z = 6, and otjv = 1-6. All 

simulations allow for a random dispersion in otjv with a a = 0.2. At fixed £, the number of photons per timestep is the same in all the 

simulations that use QSO Model II, independent of ouv- 

b The "LLS" column specifies the method for capturing the dense absorbers, as discussed in Appendix lAl 

c iso = isotropic emission 
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Fig. 3. — Volume-averaged global history of ThoU in s 1 , 
kilo-Kelvin, and ShoIHi in six of our simulations. 
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discussed in this paper, and Figure [3] plots the volume- 
averaged history of ThcH, T, and XHein m srx °f these 
simulations. The simulations in Figure [3] use different 
prescriptions to include the QSOs and to capture the 
high column-density absorption systems. Yet, with a mi- 
nor amount of tuning, the end of Hell reionization occurs 



near z ~ 3 in all of the simulations. 16 The known pop- 
ulation of quasars is able to ionize the Hell in the Uni- 
verse by 2.5 < z < 3.5, approximately the redshift range 
where several observations suggest that Hell reionization 
is ending. In the end, 2.2 photons per helium atom are 
required to reionize the Hell to XHcii = 0.95 in simula- 
tion D2 (which employed no filtering method) and 3.4 
photons per helium atom are required in simulations LI 
and L3 (which used filtering methods A and B). 

The middle and top panels in Figure[3]show the history 
of T and ThoII- In these simulations the temperature of 
the IGM peaks near the end of Hell reionization, and 
Fueii increases during this process. The jaggedness of 
rHeii(^) owes to sample variance. 17 

The remainder of this section discusses our simulations 
in detail. But first, it is important to highlight a draw- 
back of our radiative transfer method. Time steps are 
set to At = 10 Myr, such that a ray will travel a dis- 
tance c At over a timestep before being held in memory. 
This introduces a characteristic scale of c At = 3 proper 
Mpc that can appear in the simulations. In particular, 
the effect of a finite timestep is sometimes apparent in 
the THeii field, through rings that show sharp gradients 
in ThcH- These small artifacts do not affect our conclu- 
sions. 

16 Note that we tune the parameter f , which normalizes the total 
ionizing emissivity, between 0.75 and 2 in the simulations presented 
here in order for Hell reionization to end at z ~ 3. Given observed 
uncertainty in the QSO luminosity function and spectral behavior 
near the Hell Lyman-limit, such tuning is acceptable. This aids 
comparison between simulations, and it also allows the simulations 
to address better the set of observations that indicate that Hell 
reionization is nearing completion at z m 3. 

17 The jaggedness is particularly acute in simulation S2 because 
it uses the shortest quasar lifetimes, it has the fewest quasars at 
high redshifts, and its brightest quasars have the same lifetime as 
its dimmest ones. 



4.1. Homogeneous Universe 

Simulation D2 is run in the 186 Mpc box using the 
positions of the halos in the N-body code and Method 
II for populating these halos with quasars. In order to 
isolate some interesting physical effects, this simulation is 
run in a homogeneous density field (Ab(x, z) = 1) rather 
than using the density field from the N-body simulation 
as in the other cases. The left-most panels in Figure 2] 
depict the Hell ionization field at several different times. 
Initially, the Helll bubbles around each QSO are isolated 
spheres with characteristic circular radii of w 20 Mpc. 18 

The edges of the Helll regions are not sharp, with the 
region over which 0.1 < iehcIH < 0.9 spanning several 
Mpc and with a lower ionization tail extending much 
further. When the Universe is roughly a few tenths ion- 
ized, the spherical HcIII regions begin to merge and form 
larger non-spherical Helll regions. In addition, the frac- 
tionally ionized regions that extend far ahead of the ion- 
izing front overlap, creating a global ionization floor. By 
the time xhcIH > 0.8 in simulation Dl, xhoIII > 0.2 ev- 
erywhere. 

Since the recombination timescale at mean density is 
shorter than the Hubble time and comparable to the 
Hell reionization timescale, relic Helll regions recom- 
bine significantly over the course of Hell reionization. 
Relic Helll regions in which the Helll has substantially 
recombined are present in the top-right panel in Figure 
|U However, these regions persist only while the Hell 
ionizing background has ThcH < «b S e « 2 x 10~ 17 s _1 
at z ss 4, where «b is the Case B recombination coeffi- 
cient. Otherwise, ionizations will be sufficient to balance 
recombinations. Relic Helll regions that have substan- 
tially recombined exist in our simulations primarily dur- 
ing the early phases of Hell reionization (xhcIII i$ 0.5). 
The fourth column in Figure [4] displays ThcH- The con- 
dition ThcH > 2 x 10~ 17 s _1 is met in most regions in 
the IGM at times when ShoII j^ 0.5. 

The panels in the second column of Figure[4]show log 10 
of the transmission in the Hell Lya forest. For a homo- 
geneous Universe, there is little transmission in the Hell 
Lya forest, even at times near the end of Hell reioniza- 
tion. The inclusion of density inhomogeneities signifi- 
cantly enhances the transmission level (Section 14. 2 [I . 

The IGM temperature field has a different morphol- 
ogy from the xjjein field (compare the third column in 
Fig. [|]with the first). Early on in Hell reionization, the 
morphologies of both the XHciii and T fields are driven 
by isolated Helll regions. The temperature of an iso- 
lated Helll region is AT ps 7, 000 K above the average 
temperature in Hell regions, with hotter regions at the 
Helll front edge and slightly colder regions inside the 
front. Once ionized, the HcIII region begins to cool adi- 
abatically (owing to cosmic expansion). Regions that 
were ionized earlier are cooler than regions that have 
been ionized more recently. The volume-averaged tem- 
perature increase of 10, 000 K during Hell reionization is 

18 This value for Rj, may seem somewhat smaller than those pro- 
jected in Figure [2] for several reasons: (1) A random slice through 
a bubble will on average have a factor of 2 smaller radius than 
the 3-D radius; (2) active quasar bubbles are smaller than the fi- 
nal bubble size - the quantity plotted in Figure [2] (the quasars in 
Method II have long lifetimes that are comparable to the Salpeter 
time); (3) a fraction of the ionizing photons are deposited ahead of 
the front and some reionize recombining atoms; (4) f < 1 in D2. 



consistent with estimates that assume that the IGM ab- 
sorbs all photons with energies less than a few hundred 
eV during Hell reionization (Section 2.3). 

The right-most column in Figure [4] displays the cumu- 
lative amount of Hell photo-ionization heating measured 
in Kelvin. In the temperature panels it is apparent that 
the coolest regions are the first regions to be ionized, 
which we justified by the fact that these regions had more 
time to cool. However, in these photo-heating panels, we 
see that the regions that are ionized first also have the 
least heat injection from Hell photo-heating. Regions 
that are ionized at later times have absorbed more of the 
hard radiation background before being fully ionized. 

4.2. Density Fluctuations 

Figure [5] is the same as Figure |4j but features simu- 
lation Dl which includes density inhomogeneities using 
the 256 3 gridded N-body field. The small-scale structure 
of Hell reionization is changed significantly by including 
density inhomogeneities (compare the left-most panels in 
Figures [4] and [5]) . However, the large-scale morphology 
of Hell reionization is similar between the homogeneous 
and inhomogeneous cases. 

The structure of the temperature fluctuations is af- 
fected the most by density inhomogeneities compared to 
the other quantities shown in Figure[5] Adiabatic heating 
and cooling from structure formation imprints additional 
small-scale features into the temperature field in simula- 
tion, and it makes more regions cooler in Dl than in D2 
because most of the volume in the IGM is underdense 
and expanding faster than if it were in the Hubble flow. 

The second column in Figure [5] shows the total Hell 
Lya forest transmission. The transmission level is signifi- 
cantly higher when density inhomogeneities are included. 
Underdense regions are responsible for most of the trans- 
mission in the z *=s 3 Hell Lya forest. Section l5~3l shows 
that simulation D2 is consistent with measurements of 
the Hell Lya forest mean transmission. 

Figure [6] presents simulation Bl, which has the same 
run specifications as simulation Dl, but in the larger 429 
Mpc box. The larger box provides a better sample of 
the structures during Hell reionization. On these larger 
scales, QSO clustering is evident in the maps. 

4.3. The Effect of Dense Clumps 

Figure [7] compares three simulations that employ dif- 
ferent methods for including the effect of dense clumps on 
filtering the ionizing radiation. Each row features a snap- 
shot from its respective simulation with ShcIII = 0.85. 
The top panels are from simulation Dl (which does not 
employ an additional filtering method), the middle pan- 
els are from simulation LI (filtering Method A), and the 
bottom panels are from simulation L3 (Method B). All 
three simulations are run on the 256 3 grid using quasar 
Model II. Simulation Dl reaches XHciii = 0.85 at z = 3.3, 
and simulations LI and L3 approach this mark at slightly 
later times, at zk 3.15, owing to the additional filtering. 
For the most part, the XHein and T fields show similar 
patterns in the three simulations, which suggests that 
our results depend weakly on the filtering method. The 
volume-averaged temperatures in the three simulations 
are within 2000 K of each other at fixed XHeii- Further- 
more, the distribution of temperatures at mean density 
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Fig. 4. — Slices through simulation D2, which assumes a homogeneous IGM (Aj, = 1). Each slice has breadth 186 Mpc and width 3.5 Mpc 
(except for the second column, which has a width of 1 grid cell or 0.7 Mpc). From top to bottom, the panels in each row feature snapshots 
at z = 4.3,3.5,3.2 and 3 respectively, corresponding to 5jh b iii = 0.1, 0.5, 0.8 and 0.99. The first column shows the spatial distribution 
of XHelll- The second column depicts log 10 of the transmission in the Hell Lya forest. The middle column shows the temperature of 
the IGM. The fourth column plots log 10 of r[i e ii, and the right-most column shows the cumulative amount of heat deposited by Hell 
photo-ionizations. White regions represent higher values than shown in the scale and black represent lower values. 



are similar between these simulations, but with the distri- 
bution in LI extending to a few thousand Kelvin higher 
temperatures than those in the other simulations. 

The right-most column in Figure [7] illustrates how the 
filtering methods affect the m.f.p. of ionizing photons, 
where the m.f.p. is tabulated by averaging the dis- 
tance that absorbed photon travel (where this average is 
weighted by the energy it deposits). The average m.f.p. 
decreases slightly from simulation Dl to LI as expected, 
and again from LI to L3. 

It is surprising that filtering from dense systems does 
not have a more significant impact, particularly on the 
temperature fields, which are sensitive to the hardness 
of the ionizing radiation. If the bubble sizes were much 
smaller than the photon m.f.p., then it would not be sur- 
prising that filtering by dense regions is unimportant - 
radiation would be filtered minimally by dense clumps 
prior to reaching diffuse neutral gas. However, the bub- 
ble sizes are comparable to the m.f.p. (Appendix A.), 
and we know that about 1/3 of the photons are absorbed 
in dense systems in LI and L3 since Hell reionization in 
these simulations requires 3.4 photons per Hell ion as 



opposed to 2.2 in simulation Dl. If we assume that the 
filtering removes the bottom third of the Hell-ionizing 
photons in energy and that the IGM is ionized by the 
remaining, harder radiation, the temperature increase of 
the IGM would be a factor of 2.4 larger for auv = 1-6 
than if the filtering had not occurred (at least if AT is 
computed in the optically thin limit; eqn. [6]). 

Why is filtering in our simulations different from this 
toy scenario? The answer is that not just the lowest 
energy photons are filtered. The effective optical depth of 
dense systems s cales as ~ -E" 1 ' 5 if the Audi distribution 
scales as —1.5 l|Zuo fc Phinnevlfi993h . This absorption 
contrasts with the absorption from diffuse neutral gas, 
for which the effective optical depth will scale closer to 
~ E~ 3 . Therefore, the spectral shape of the ionizing 
radiation is not as altered by filtering by dense systems as 
by diffuse IGM gas, and it is certainly not as hardened by 
dense systems as in the above toy model. Furthermore, 
equation (|7|) shows that the amount of heating depends 
weakly on the spectral index of the incident spectrum. 

Much of the filtering that hardens the radiation field 
in the Hell photo-heating panel in Figure [7] occurs in 
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Fig. 5. — Same as Fig. |4]but using simulation Dl (which includes density fluctuations). The panels in each row are chosen to match the 
^Helll fractions quoted in the caption of Fig. [4] 



the diffuse gas rather than in dense clumps. This is why 
the temperatures in simulation D2, for which there are 
no density fluctuations, are similar to the other simula- 
tions. Even in the Helll bubbles in Figure [3 there are 
significant amounts of neutral gas which do not reside 
in very overdense locations. Furthermore, photons that 
penetrate into neutral regions will typically encounter a 
higher total Hell column density in diffuse gas than in 
dense systems. These photons are responsible for heat- 
ing the hottest regions (which are the last regions to be 
ionized). 

Figure [5] compares the globally averaged spectrum in 
simulations Dl (no filtering) and LI (filtering Method 
A). This global average is computed by tabulating the 
incident spectrum in randomly selected cells in the simu- 
lation volume and then averaging. (The global spectrum 
in L3, which is not shown, is similar to this spectrum in 
LI.) Reionization ends (xhcIU = 0.95) at z = 3.3 in Dl 
and at z = 3.1 in LI. The average spectrum has a similar 
shape, particularly at the beginning of Hell reionization 
when the impact of filtering is minimal. At later times, 
the difference is a factor of < 2, and the slope is not 
significantly changed. Near the end of Hell reionization, 
the amplitude of the spectrum evolves by order unity in 
both simulation Dl and LI on timescales of Az s=a 0.2. 



This fast evolution suggests that diffuse gas is playing a 
role in limiting the mean free path in these simulations 
even at the end of Hell reionization. 

After Hell reionization, the filtering method affects the 
spectrum more because the m.f.p. is limited by dense 
systems and not by diffuse gas. The value of ThcH, which 
in a large part determines the transmission in the Hell 
Lya forest, is proportional to A m f p . The photo-ionization 
rate at the end of the simulation is 2 — 4 x 10~ 15 s _1 in 
LI and L2 whereas it is 7 x 10~ 15 s _1 and increasing 
steadily when simulation Dl was terminated at z = 2.9. 

4.4. Soft X-ray Photons 

Simulation XI is the same as LI except that photons 
with m.f.p. larger than 3.5 times box size, background 
photons with roughly E n > 500 eV, are not included in 
the calculation. How much heating and ionizations do 
these high-energy photons contribute? When the Uni- 
verse is 90% ionized in Hell, model XI is 1400 K cooler 
than model LI. Furthermore, the fraction of ionizations 
that arise from photons with E 1 > 500 eV is small. The 
heating from the background is largest in the regions that 
are ionized last, the regions that are exposed to the back- 
ground for the longest period of time. The background 
heating adds an additional several thousand Kelvin in 
these regions. 
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Fig. 6. — Slices through simulation Bl, which employs the same 
run specifications as simulation Dl (Fig. QJ}, but uses the 429 Mpc 
box. Each panel shows a 429 X 429 X 8.4 Mpc slice through this 
simulation volume. From top to bottom, xjiolll = 0.1, 0.5, and 0.9. 

The relative unimportance of heating from photons 
with E~ t > 500 eV is reassuring, because we do not 
model the obscured quasar contribution to the ioniz- 
ing background. Models of the obscured contribution to 
the diffuse X-ray background find that obscured quasars 
(defined as having QSOs with obscuring columns of 
JVhi > 10 21 cm~ 2 ) begin to do minate the backg round 
above approximately 2 keV fe.g.. lGilli et alj |2007). 

4.5. QSO Model 

The major uncertainties associated with modeling 
QSOs are their lifetimes, spectra, and whether their emis- 
sion is beamed. Figure [9J considers 3 QSO models, to 
illustrate how these uncertainties influence our results. 
From top to bottom, the panels represent snapshots with 
xneiii ~ 0.1, 0.75, and 0.9. Beamed emission (see simula- 
tion SI in Fig. [9]which uses the beaming model described 
in Appendix E) enhances the structure of the ionization 
field earlier on in Hell rcionization. Similarly a light bulb 
model with tqso = 10 Myr (simulation S2 in Fig. [§]) re- 
sults in more relic Helll regions compared to our fiducial 
case and in slightly smaller bubbles. This fact owes to 
the shorter quasar lifetime compared to the fiducial case. 



The IGM temperatures in S4 and S4b, which have a 
harder UV spectral index of 0.6 compared to the fiducial 
value of 1.6 (cho sen to better agree with the results of 
IScott et al.ll2004l ), are significantly higher than the other 
simulations analyz ed in this paper. T his result agrees 
with that of iTittlev k. Meiksinl (|2007h . who also find a 
strong dependence of the temperature on the spectral 
index of the sources. In S4 (right panels in Fig. [9]), 
the hottest regions are typically the least ionized. In 
the other simulations featured in Figure [H the hottest 
regions are typically the most ionized until XHein ^ 0.8. 

Why is the morphology of the temperature so differ- 
ent in S4 whereas the aJHeiii field is similar to the other 
simulations (except for a higher ionization floor)? For 
<5uv = 0.6, the average energy photon that is absorbed 
is much higher than in the fiducial model, injecting more 
heat. However, near quasars it is still the softest photons 
that are absorbed, resulting in similar structures in the 
iCHeiii field compared to the fiducial model. The hardest 
photons (which are more plentiful in S4) are absorbed in 
the neutral regions outside of the Helll bubbles, leading 
to the inverted behavior between XHeiii an( l T. 



4.6. Power Spectrum 

We have seen that 50 Mpc fluctuations in the temper- 
ature and XHoiii are present during and after Hell reion- 
ization. These fluctuations modulate the level of absorp- 
tion as well as the amount of small-scale power in the HI 
Lya forest. Here, we use power spectrum statistics to 
quantify the properties of these spatial fluctuations. 

The curves in Figure [TUJ are calculated from simula- 
tion Bl, which uses the 430 Mpc box, and represent the 
power spectrum of ^Heiii (top panel), ST = T/T—l (mid- 
dle panel) , and the cross correlation coefficient of rEHeiii 
and 5b (bottom panel) from simulation Dl, which is de- 
fined as r = Pg s h /\Ps s .fk«J -1 i where 
S XlIoII1 = xhoIii/ShcIii ~ 1- We have checked that the 
power spectrum in Bl agrees well with that in Dl, which 
is identical to Bl except run in the 429 Mpc box. The 
only scales where the two simulations do not agree corre- 
spond to k > 2 h Mpc~ , where there is an upturn owing 
to shot noise in Bl that is not present on these scales in 
Dl. 

The XHoiii fluctuations peak on ~ 50 Mpc scales, even 
early in the Hell reionization process when XHeii sw 0.1, 
with the peak in the amplitude of the fluctuations oc- 
curring at JCHein ~ 0.5. This peak scale is remark- 
ably constant throughout Hell rcionization. This con- 
trasts with the evolution of the power spectrum in most 
models of hydrogen reionization, in which the growth of 
large HII regions around clusters of sources leads to the 
power shifting to larg er scales with increasing xm (e.g., 
iMcQuinn et al.ll2007h . 

The fluctuations in the temperature field are shown 
in the middle panel in Figure [TU] at several xhcIU- F° r 
comparison, the thin blue dotted curve is from a simula- 
tion where Hell reionization does not occur at z < 6 and 
that is initialized to have 7 — 1 « and T = 10 4 K at 
Z = 6. This curve is at the same redshift, z = 4.0, as the 
^Heiii = 0.5 curve. This illustrates that Hell reionization 
produces temperature fluctuations on much larger scales 
than structure formation. 

At k < 0.3 Mpc -1 , the fluctuations from patchy Hell 
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Fig. 7. — Three rows use snapshots at XHelll = 0.85. The top panels are from a simulation that does not employ an additional filtering 
method, the middle panels are from a simulation that uses filtering method A, and the bottom panels are from a simulation that uses 
filtering Method B. 
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Fig. 8. — Globally averaged radiation spectrum at various times 
during our simulations, computed by tabulating the incident spec- 
trum in many randomly selected cells and then averaging. The 
thin curves are for simulation Dl and the thick curves are for LI. 
LI differs from Df in that it uses filtering method A. 



rcionization are the dominant source of the T fluctua- 
tions. The scale of the peak in these fluctuations is 
comparable to the Hell bubble sizes early on. How- 
ever, when SHoiii = 0.3, 0.5 and 0.97, large fluctuations 
are present on all scales captured in our box. Interest- 
ingly, the amplitude of large-scale temperature fluctua- 
tions are largest at the beginning of Hell rcionization, 
with (AT/T) rms = [fc 3 AT/(27r 2 T)] 1 / 2 = 0.2. Fluctua- 
tions in tem peratu re are observable in t he HI Lya forest 
(|Zaldarriagall2002l; iTheuns et alj|2002al) . 

The degree of correlation between the ionization and 
density fields reveals how important QSO clustering is 
in shaping the morphology of Hell reionization. The 



bubbles begin to correlate with the density field on 50 
Mpc scales (bottom panel in Fig. [T0|) . On smaller 
scales, ccHein is anti-correlated with 6b- The level of anti- 
correlation increases with xhoIH as the Helll regions en- 
compass more dense neutral systems. 

5. OBSERVATIONAL IMPLICATIONS 

This section discusses how Hell reionization affects 
observations of the IGM. The observables that are ad- 
dressed are the IGM temperature, the T-A?, relation, and 
the mean transmission in the HI and Hell Lya forests, 
all of which can be estimated from high-redshift quasar 
spectra. 

5.1. Temperature Evolution 

Figure [TT] plots the evolution of the temperature PDF 
in simulations D2, LI, and S3, tabulated from grid cells 
with —0.1 < 6b < 0.1. Simulation D2 has no density 
fluctuations and auv = 1-6; simulation LI includes den- 
sity fluctuations, auv = 1-6, and filtering Method A; and 
simulation S3 is the same as LI but with auv = L2. 

The evolution of the PDF in simulation LI is charac- 
teristic of our simulations with auv = 1-6. At times in 
LI for which x Hc iii < 1, most of the IGM has T « 1 x 10 4 
K, except ionized regions which have T s=s 1.6 x 10 4 
K. The result is a bimodal temperature PDF (see the 
^Heiii = 0.1 curve in the middle panel in Fig. Ill[) . 
As Hell reionization proceeds, the PDF shifts to higher 
temperatures. If only the insides of Helll regions were 
heated, this PDF would become broader rather than 
move towards higher T in the manner seen in Figure 
111! implying that hard photons are streaming far from 
their sources prior to being absorbed. The bimodal- 
ity of the PDF decreases with increasing ShcIII until 
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bers to iMcDonald et al.l (|2001l ) from measurements of 
the HI Lya for es t pow er sp ectrum. Although , both 
IMcDonald et al.1 (|2001l ) and iZaldarriaga et alj (J2001D 
measured the temperature in three coarse bins centered 
at z = 2.4, 3.0 and 3.9, possibly obscuring a sudden tem- 
perature change. 
The bottom panel of Figure [T^] plots the Schav e~et al.l 



Fig. 9. — Morphology of Hell reionization for three possible source models. In SI the emission from the quasars is beamed as described 
in the text, in simulation S2 the quasars have a light bulb behavior with tqso = 10 Myr, and in simulation S4 the quasars have a hard 
spectrum with <5tjv = "-6- From top to bottom, the panels represent snapshots with aiHelll ~ 0.1,0.75, and 0.9. Note that temperature 
panels for simulation S4 have an adjusted range that extends to 40,000 K compared to 25,000 K for the other simulations. 

ShoIH = 0.6, at which point the PDF is fairly Gaussian. 
For Snein > 0.6, the PDF again becomes skewed as the 
harder photon background builds up and starts to appre- 
ciably heat regions far from quasars. The temperature 
PDF at the end of Hell reionization is broad, extending 
between 12,000 K and 35,000 K in simulation LI. The 
temperature of a gas parcel at the end of Hell reioniza- 
tion is a result of many factors, including the redshift(s) 
it was ionized and the hardness of the radiation field that 
ionized the parcel. 

The evolution of the temperature PDF in simulations 
D2 and S3 are similar to that in LI. The PDF in D2 
does not extend to as high temperatures as in LI. This 
implies that density inhomogeneities contribute to the 
temperature of the hottest gas parcels. The evolution of 
the PDF is also similar in S3 to LI, except the harder 
spectrum in S3 results in higher temperatures. 

iSchave et all (|2000 ) and Ricotti et all ((2000) estimated 
the temperature of the IGM from the widths of the 
narrowest lines in the HI Lya forest, calibrating their 
measurements with numerical simulations that essen- 
tially assumed a power-law T-Ab relation at relevant 
Ab. These studies claimed to have detected a sud- 
den increase in the IGM temperature between z rj 3.5 



and 3. A subsequent ana l ysis u sing a similar method- 
ology by IMcDonald et al.1 {2001) did not detect a sud- 
den increase in temperature, but rather a temperature 
at mean density of T n » 1 7,000 ± 2000 K at z = 2,3, 
and 4. IZaldarriaga et all (|2001h derived similar num- 



(|2000h . lRicotti et al.l (|2000f ). and lMcDonald et al.1 (2001) 
measurement values for To as well as the evolution of this 
quantity in simulations LI, Dl, and S4b (thick curves). 
In general, our calculations over-predict the measured 
values of Tq. However, since these observations look for 
the narrowest lines in the forest, it is probable that they 
are most sensitive to t he coolest temperatures at a given 
density as is argued in iFurlanetto &; Oh! (|2007l ) . A simi- 
lar argument applies for measurements using the HI Lya 
forest power spectrum. The thin curves bracketing the 
thick ones are ±1 s.d. of the mean To. Note that the 
agreement with measurements is better if we compare 
with the —1 s.d. curves rather than the mean Tq curve. 
A more detailed study is required in order to determine 
whether these observations are consistent with our sim- 
ulations. 

Figure [T^] also includes two curves for the case in which 
Hell reionization happens at z > 6 and the simula- 
tion is initialized with temperatures of T(Af,) = 10 4 
K and T(A b ) = 2 x 10 4 K at z = 6 (thin cyan dot- 
dashed curves). By z « 3 both of these curves have 
asymptoted to a temperature of rs 7, 000 K as expected 
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Fig. 10. — Power spectrum of xjielll (top panel), ST (middle 
panel), and the cross correlation coefficient of 3>H e lll arl d <5fc (bot- 
tom panel) in simulation Bl. The thin blue dotted curve in the 
middle panel is from a simulation without Hell reionization. 
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Fig. 11. — PDF of the temperature tabulated from grid cells 
with —0.1 < <5(, < 0.1 in simulations D2, LI, and S3. 




o 

M 





30 
25 
20 
15 
10 



LI 

Lib 



-^H- 



+ 



X McDonald et al, (2001) - 
.. V, Ricotti et al. (2000) 
'**•♦.. Schaye et al. (2000) 
Dl 
LI 
V \ S4b 



+ 











2.5 



3.5 



4.5 



5.5 



Fig. 12. — Top Panel: Best-fit power-law index of the T-A b 
relation in two simulations (thick curves), as well as its evolution 
in simulations without az~3 Hell reionization (thin curves). The 
markers with error bars are measurements of this quantity in the 
literature, and the references for these measurements are given in 
the key in the bottom panel. Bottom Panel: Evolution of best- 
fit To in three simulations (thick curves) as well as ±1 s.d. in 
this quantity (the thinner curves with the same line style). The 
thin cyan dot-dashed curves that are increasing with z represent 
simulations in which Hell reionization occurs at z > 6. Also shown 
are measurements of To . 



(|Hui fc Gn cdin 1997J). The curves that represent the sim- 
ulations with Hell reionization (at least in models Dl and 
LI) are more consistent with the measurement points 
than these curves. 

5.2. Temperature- Density Relation 

The temperature of the gas depends primarily on the 
competition between photo-heating and adiabatic cool- 
ing. These processes lead to a power-law relationship 
between T and A& for unshocke d gas at times suff i- 
ciently after a reionization epoch (|Hui fc Gnedinlfl997|) . 
This power- law index asymptote s in time to the value 
7 - 1 = 0.6 dHui fcCxnedinl [1991 . 19 

However, Hell reionization can lead to a more compli- 
cated relation between T and Af, than a single power-law. 
Just after an instantaneous and homogeneous Hell reion- 
ization process, the equation of state 7 — 1 would equal 
zero (an isothermal IGM), at least if the temperature af- 
terward is much larger than before. However, the finite 
duration of Hell reionization causes 7 — 1 to deviate from 
zero, and the duration as well as the inhomogeneity of 
the photo-heating introduces scatter into the T-Af, re- 
lationship. Regions that are ionized earliest have more 
time to cool, and regions ionized by the most filtered 
radiation are the hottest. 

Figure [13] shows the evolution of the average T-Ab re- 
lation in simulation Dl, which uses Quasar Model II, at 
times for which ShcHI = 0, 0.5, 08, and 1 (thick solid 
curves) . The slope of this relation (which was initialized 
with 7 — 1 = 0.3 at z = 6) does not change substantially 
throughout Hell reionization. While the value of 7— 1 ap- 

19 This value for 7 — 1 is just slightly smaller than if only adia- 
batic cooling is included, which would yield 7 — 1 = 2/3. 
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Fig. 13. — T-Aj, relation in simulation Dl. This simulation is 
initialized with 7 — 1 = 0.3 at z = 6.1. The thick solid curves 
are the mean T(A;,), and the contours enclose (33, 67, 90, 99, 99.9) 
percent of the grid cells — encompassing the regions that have the 
highest densities in logT - log A;,. 



pears fairly constant, the dispersion in this relationship 
grows during Hell reionization. The contours enclose 
(33, 67, 90, 99, 99.9) percent of the cells in the simulation 
grid - encompassing the regions that have the highest 
densities in logT - log At . 

The top panel in Figure [T2l shows the measured value 
of 7 — 1 at various times in simulations LI and Lib, which 
use QSO Model II and filtering Method A. 20 The only 
difference between these two simulations is that LI is ini- 
tialized with 7 — 1 = 0.3 at z = 6 and Lib with 7 — 1 = 0. 
In both of these cases, as well as in all of the simulations 
listed in Table [TJ the T-Aj, relation does not become in- 
verted or isothermal (i.e., 7 — 1 < 0). This result owes 
to the large QSO bubbles being essentially uncorrelated 
with Aj,; no value of Aj, is ionized preferentially at a 
given time during Hell reionization. Unless some A& is 
ionized preferentially, then the inequality 7 — 1 > must 
hold. Our simulations yield 7 — 1 m 0.3 during the bulk 
of Hell reionization. After Hell reionization, this index 
steepens. Interestingly, the evolution of 7(2) — 1 in sim- 
ulations S4, which has auv = 0-6, is nearly identical to 
this in LI, even though the IGM in S4 is heated to a 
much higher temperature. 

Also shown in Figure [12] are measurements of 7 — 1 
in the literature. The different estimates are not con- 
sistent with one another, and generally have large error 
bars. These measurements are calibrated using hydrody- 
namic or hydro-PM simulations that do not account for 
inhomogeneous heating owing to Hell reionization. 21 

Previous studies have found differing results concern- 

20 We fit a power-law to the mean value of T(A(,) in bins span- 
ning 0.5 < A(, < 2 to obtain 7— 1. We find that the same procedure 
applied to the median of T(A(,) rather than the mean or to a more 
restricted range in A;, yields similar results. 

21 Also, note that the Schavc et al. (2000) data points - which 
are the most discrepant with respect to our predictions - are 
each taken from half of the Lya forest from a single quasar spec- 
trum (i.e., approximately a 250 Mpc skewer). Because of the 
large structu res during Hell reionization, it is conceivable that the 

ISchave et al.l l)2000fl estimates do not sample enough regions to be 
representative. 



ing the effect of Hell reionization on the T-Afc rela- 
tion. As in our simul a tions , the T-A& relation in the 
studies of lGleser et"aH (|2005| ) and iPaschos et all (|2007l ) 
is not inverted. I n fact, the predictions for 7 — 1 in 
iGleser eta l. (2005) are qui t e steep with 7 — 1 > 0.38. In 
contrast, iFurlanetto fc Obi (|2007t i found that the T-A b 
relation could become inverted by Hell reionization and 
that the mean relation could deviate significantly from a 
power-law in their excursion set based model. This result 
owed to the underdense regions being ionized later dur- 
ing Hell reionization and, therefore, being the hot test. 
In a toy case considered in IFurlanetto fc Obi (|2007l ). in 
which there are no correlations between the QSO bubble 
and the density field, the evolution of 7 — 1 during Hell 
reionization is similar to what we find. 

The no-corr el ation approximation in the 
IFurlanetto fc Ob] (|2007l ) toy model is reasonable 
during Hell reionization by quasars. The magnitude 
of the correlation between Helll bubbles and fluc- 
tuations on the Jeans scale is relevant for how Hell 
reionization affects the T-Ab relation. The linear-theory 
correlation coefficient between the overdensity in a 
region of 10 12 M that hosts a QSO and a Jeans mass 
region separated from it by 20 Mpc (where we take 
Mj = 10 10 M Q ) is r = 0.02 (r = 0.008 for a 30 Mpc 
separation). The value of r would be even smaller if 
non-linearities are taken into account. Therefore, the 
probability distribution function of 8b smoothed on 
the Jeans scale within the Helll region of a QSO is a 
Gaussian (in linear theory) centered at (5b) = rN a , 
where N a is the quasar overdensity in units of the s.d. 
in 8b- The variance of this Gaussian is reduced by 1 — r 2 
from the volume-averaged variance. Therefore, for the 
small r values quoted above and reasonable values of 
N a , a Jeans mass clump in a quasar Helll region has 
a nearly equal chance to be overdense as it does to be 
u nderdense. 

iBolton et al.1 (|2008al ) and other studies have postulated 
that radiative transfer effects can result in an inverted T- 
Ab relation. The idea is that the radiation which makes 
it out into the underdense voids is highly filtered, heat- 
ing these regions to higher temperatures. However, this 
picture does not apply when quasars are the sources that 
ionize the Hell. Voids on the Jeans scale, the scale that 
is most relevant for the Lya forest transmission, are not 
strongly anti-correlated with the giant Helll bubbles. In- 
stead, radiation that is filtered as it traverses a QSO bub- 
ble has an approximately equal chance of ionizing a void 
or filament at the bubble edge. Therefore, the filtering 
leads to scatter in the relation between T and A&, but 
not to an inverted relation. If the Hell-ionizing sources 
were more numerous and much dimmer than QSOs, in 
this case an in verted relation is poss ible. 

Interestingly, IBolton et al.l (|2008af ) found that the flux 
probability distribution from the z ~ 3 HI Lya forest 
favors an inverted T-A& relation to accommodate the 
measured number of high transmission pixels. We have 
shown that the observed population of QSOs cannot be 
responsible for an inverted relation. However, we think 
that Hell rei o nizatio n by QSOs c ould still reconcile the 
IBolton et all (|2008af) result. The IBolton et all (|2008af) 
analysis did not consider physically motivated models 
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Fig. 14. — Evolution of the effective Hell optical depth. The 
error bars on most of the simulation curves bracket the values in 
which thcII eff falls 90% of the time when measured over a skewer 
of length 190 Mpc. The measurement points are described in the 
text. 

for the dispersion in this relation. 22 The large dispersion 
in T(Ab) that we find will act to increase the number of 
una bsorbed pixe l s and could potentially accommodate 
the iBolton et a l. (2008a|) result without the need for an 
inverted T-A;, relation. 

5.3. Mean Hell Lya Forest Transmission 

There are currently only a handful of Hell Lya for- 
est sight-lines. The number of such sight-lines is limited 
because a bright quasar is required to achieve adequate 
sensitivity with existing instruments and because a sin- 
gle HI Lyman-limit system can absorb much of the flux 
at wavelengths where Hell Lya absorption occurs. The 
number of Hell Lya forest sight-lines is expected to in- 
crease once the Cosmic Origins Spectrograph is installed 
on Hubble, which is currently scheduled for spring of 
2009. 

In this section, we calculate T e s t n c u along skewers 
through our simulation box. These spectra, as well as 
those discussed in Section 15.41 are generated by grid- 
ding the N-body particles on a 256 3 mesh in the 186 
Mpc box and using the temperature and ionization fields 
from the simulations. Our predictions for the HcII-Lya 
transmission are reliable to the extent that our dark mat- 
ter simulations capture the gas density PDF (especially 
for 5b < 0, where most of the transmission occurs). In 
Appendix [DJ we show that the dark matter PDF mea- 
sured from our grid reasonably reproduces the expected 
gas PDF. Our predictions for the mean Hell Lya trans- 
mission change by only < 20% when we calculate this 
statistic from the 512 3 gridded N-body field rather than 
the 256 3 one. This agreement is reassuring and suggests 
that our predictions are robust. 

Figure Q3] plots Tcf^Hcii from our simulations and 
from the most constraining observations, where Tef^Hen 

22 They did consider a model with 7 — 1 m 0.3 in which this 
scatter was increased artificially by adding a Gaussian dispersion 
to the particle photo-heating rates, and they found that this model 
resulted in essentially no c hange in the flux PDF over a model with- 
out scatter. However, the Bolton et al. (2008a) model with scatter 
resulted in similar 90% contours in the T-Aj, plane compared to 
their mo del without scatter (see their Fig. 3 and compare with our 
Fig. 113ft . The width of their 90% contours in the T direction are 
factors of several smaller than what we find. 



ExD HS1157+3143 llReirners et all l2005h PKS1935 
692 (jAnderson et al.l |l999|), Q0302-003 (lHejffl _et_aL| 
120001 and SDSSJ2346-0016 (jZheng et al.ll2004h . The 
vertical error bars on these measurements are just sta- 
tistical, and the horizontal error bars signify the redshift 
range over which each of the values were estimated. 23 

The curves in Figure [H] are Tef^Heii from four simu- 
lations, for which each marker value is estimated from 
1000, 186 Mpc random skewers through the simulation 
volume. The error bars on most of the simulation curves 
bracket the values in which Teff^eii falls 90% of the time 
when measured along a single skewer. Note that the 
measurement points were taken from skewers that are 
comparable in length to the 186 Mpc box. There should 
be an additional cosmic variance error on these points 
that is similar in magnitude to the error bars on the sim- 
ulation curves. 

Our simulated predictions for Tef^Heii^) provide mod- 
erate agreement with the measured values. The evolution 
of Tof^Hoii^) in simulations Dl and Lie is more consis- 
tent with three of the four measured points at z < 3 than 
the T e s,Koii{z) in simulation LI. None of the simulations 
are consistent with the z = 2.7 measurement (although, 
simulation Dl was terminated at z — 2.9 and would have 
been the most consistent). 

The differences in Teff^Hcii^) among the simulations 
stem in part from reionization ending at z — 3.15 in 
Dl and z = 3.45 in Lie, whereas it is completed at 
a later time, at z = 2.95, in simulations LI and L3. 
(The end of Hell reionization is defined here as when 
SHein(zcnd) = 0.95.) Another difference is that, at fixed 
ShcIH, the value of rueii is smaller in simulations LI, 
Lie, and L3 than it is in Dl because the former sup- 
plement the resolution by using filtering Method A or B 
while the latter does not employ additional filtering. A 
larger value for Fueii results in more transmission. 

Simulations LI and L3 differ only in the filtering 
method that is used, and this also leads to differences in 
the predictions for T ff : Hcii(z) primarily after Hell reion- 
ization. The uncertainties in the filtering of absorption 
systems can be remedied with improvements to the radia- 
tive transfer algorithm and by using high resolution gas 
simulations. Even though some theoretical uncertainty 
remains, the different trends among the simulations sug- 
gest that T e ff,Hcii(z) could be useful to distinguish be- 
tween different Hell reionization models. 



5.4. Mean HI Lya Forest Transmission 

By increasing the temperature of the IGM, Hell reion- 
ization changes the photo-ionization state of the inter- 
galactic hydrogen (because xhi ~ a a n e /Tm, where n e 
is the electron number density and a\ ~ T~ - 7 ), thereby 
affecting the tran smission pro perties of the HI Lya for- 
est. Interestingly, Bernardi et al.l (|2003) found evidence 
for a depression in the mean transmission at z ~ 3.2 
with width Az « 0.3. A similar depression was later 

23 Often the span over which a measurement of T c g Hell i s done 
is chosen such that it encompasses a very dark (or bright) region 
in the Hell Lyce forest. Unfortunately, this practice significantly 
biases the observational estimates. 
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Fig. 15. — Evolution of the effective HI Lyo forest optical depth. 
The red solid points with error bars are the measurement of r e fj hi 
from Fauchcr-Gigucrc et al. (2008a), and the black crosses are the 
estimates of Bcrnardi ct al. (2003). The features that arc present 
at z m 3.2 in the [Bcrnardi ct al. (2003) and Fauchcr-Gieucrc ct al. 
(2008a) data cannot be explained by Hell reionization in our simu- 
lations. Any disagreement between the data points and the model 
curves does not imply that the model is disfavored, and instead it 
may indicate that an incorrect form for rjji(z) was assumed. 



discovered bv IFaucher-Giguere et all (|2008a). 24 

The red points with s o lid erro r bars in Figure fT5l are the 
IFaucher-Giguere et al.l (|2008af) measurement of r e fr,Hi, 
which is defined as minus the log of the HI Lya for- 
est transmission. Th e black crosses with errors are the 
iBernardi et al.l (|2003f ) measurement. All errors bracket 
the 68% confidence level region. Neighboring points in 
the IBernardi et al.l (|2003h measurement are correlated. 
The dip in absorption in these measurements that has 
been the subject of attention is at z « 3.2 and has width 
Az < 0.4. A dip i s pres ent with higher significance in 
the IBernardi et all (|2003 ) measurement. 

The theoretical curves are calculated from our simula- 
tions in the manner described in Section 15.31 except for 
HI Lya absorption. These calculations assume a uni- 
form HI- ionizing background with Tni(z) — 10 -12 s _1 . 
The effective optical depth in simulation LI is represen- 
tative of most of our runs, and Tor^Hi in simulation S4b 
is smaller than is typical because of the higher IGM tem- 
peratures achieved in this simulation. The dashed curves 
give the evolution of the transmission for a simulation 
that is initialized with 7 — 1 = and To = 10 4 K or 
2 x 10 4 K at z = 6 and in which Hell reionization does 
not happen at z < 6. Any disagreement between the 
data points and the model curves does not imply that 
the model is disfavored. The form of Tm(z) can be ad- 
justed to yield any form for TeffHlO^)- 

Our predictions for Teff,Hi(.z) show no feature owing to 
Hell reionization. Hell reionization takes place over a 
much longer interval in the simulations than the Az < 
0.4 of the features seen in these measurements. The dura- 
tion of Hell reionization in our simulations is determined 
primarily by the density of L ~ L. t quasars, which over 
the relevant redshift interval is well constrained by cur- 
rent observations. In addition, even once Hell reioniza- 



24 The Fauchcr-Giguerc ct al. (2008a) estimate uses 84 high res- 
olution Keck and Magellan spectra, whereas the Bcrnardi ct al. 
(2003) measurement is from 1061 lower signal to noise and lower 
resolution Sloan spectra. 



tion is complete, it takes Az ~ 1.5 after Hell reionization 
for the gas at mean density to cool by a factor of two. 
Therefore, cooling after Hell reionization is unlikely to 
p roduce any sharp feat ure in T er r,Hi- 

iTheuns et al.l (|2002bf ) found using hydrodynamic simu- 
lations that T ffjji can evolve significantly over a shorter 
interval than the time over which the gas cools. They 
explained this effect as owing to velocity gradients that 
are induced by the heating. They claimed that these 
gradients broaden absorption lines, creating the rise in 
Teff.Hi- Furthermore, they found that that the rise after 
Hell reionization in their simulatio ns is similar to the rise 
seen at the low redshift end of the IBernardi et all (|2003f ) 
feature. This effect is not captured by our method, which 
misses the backreactio n of the photo-hea ting on the gas 
density. However. I Theuns et al.l (| 2002bl ) attributed the 
decrea se in r e ff,Hi on the other side of the IBernardi et all 
(2003) feature to heating from Hell reionization, which 
they modeled as a quick, homogeneous process. In our 
simulations, Hell reionizatio n takes place over Az > 1 
compared to Az ~ 0.1 in the lTheuns et all (|2002bf ) sim- 
ulations, and, therefore, any trends in r e ff,Hl owing to 
heating must be broader in our picture. 

6. CONCLUSIONS 

We have run a set of simulations of Hell reionization 
to understand the structure of Hell reionization and its 
effect on several observables. We find that for a late 
reionization of Hell by quasars: 

• The popular assumptions that the HcIII ionization 
fronts are sharp and that there is uniform heat- 
ing within the front (and no heating outside of it) 
do not yield a realistic model for Hell reionization. 
While the ionization fronts are still fairly localized, 
hard photons stream far from their sources and are 
absorbed ahead of the front. These photons inject 
at large distances a significant fraction of the en- 
ergy radiated above the Hell Lyman- limit. 

• The average temperature at the mean density is 
increased by w 12, 000 K over the average temper- 
ature of the gas in the absence of Hell reionization 
for ctuv = 1-6. This temperature increase is con- 
sistent with estimates that assume that the IGM 
absorbs all photons with energies less than a few 
hundred eV during Hell reionization. Regions that 
are ionized last are ionized by the hardest radia- 
tion, reaching T > 30, 000 K in our fiducial model. 
If the spectral index of the QSOs is different from 
our fiducial value of quv = 1-6, the average tem- 
perature can be be significantly altered. 

• Poisson fluctuations in the number of QSOs rather 
than their spatial clustering shape the structure 
of the ionization and temperature fluctuations on 
< 50 Mpc scales because rare L ~ T* quasars ionize 
the Hell. Hell reionization produces large tempera- 
ture fluctuations on 50 Mpc scales, and it results in 
the ionization and Hell photo-heating fluctuations 
being essentially uncorrelated with the density fluc- 
tuations on the Jeans scale, the scale most relevant 
for HI and Hell Lya forest statistics. 

• Measurements of the z ~ 3 forest suggest To ~ 
20, 000 K. Without invoking an exotic heating 
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mechanism, a late Hell reionization epoch is re- 
quired to produce this temperature. The amount 
of additional heat provided by Hell reionization in 
our simulations is enough to produce the inferred 
temperatures. 

• To the extent that there is a T-Ab relation, Hell 
reionization by quasars leads to a temperature- 
density equation of state of 7 — 1 -0.3 for 0.1 < 
^Heiii ^ 0-9. Hell reionization by QSOs cannot re- 
sult in an inverted relation (7 — 1 < 0) as has been 
claimed. 

• Our simulations of a z ~ 3 Hell reionization pro- 
cess produce a similar evolution in the Hell Lya 
mean transmission to what has been estimated. 
Better observations of the mean transmission and 
its scatter will be able to rule out models presented 
here. 

• In our simulations, the heating from Hell reioniza- 
tion by quasars is unable to produce any semblance 
of the observed depression in the z « 3.2 HI Lya 
forest opacity. The Hell reionization process is too 
extended in redshift (Az > 1) to be responsible for 
this narrow feature (Az < 0.4). 

The morphology of Hell reionization is consider- 
ably different than th at of hydrogen reionization (e.g., 
iMcQuinn et alJ 120071 ). During hydrogen reionization, 
current models find that the growth of ionized bubbles 
is more collective; hundreds or even thousands of galax- 
ies within a bubble contribute to its growth, leading to 
the bubble structure tracing the distr ibution of galaxie s 
and to tens of Mpc HII regions (e.g., IZahn et al.ll2007|) . 
During Hell reionization, the growth is more stochastic. 
Regions that happen to host a "nearby" quasar (within 
~ 30 Mpc) are ionized by that quasar. Because the QSO 
bubbles are so extended, the ionized structures are larger 
than during HI reionization. 

Another significant difference stems from the m.f.p. of 
the ionizing photons to be absorbed in diffuse gas. The 



spectrum from quasars is harder than that of stars - our 
best guess for what ionizes the hydrogen - the number 
density of helium is 70 times smaller at z — 3 than hy- 
drogen at z = 6, and the cross section is 4 times smaller. 
These three factors conspire to make the typical m.f.p. 
for a Hell ionizing photon megaparsecs rather than kilo- 
parsecs, as it is for HI ionizing photons during hydrogen 
reionization. We have seen that some Hell ionizing pho- 
tons free stream far from their sources, partially ionizing 
and heating up these regions. If stars reionize the hydro- 
gen, ionization and heating occurs within HII bubbles. 

A definitive identification of the redshifts of Hell reion- 
ization will place constraints on the sources that produce 
> 4 Ry photons and will aid studies of the HI Lya forest. 
The data from previous observations, if analyzed prop- 
erly, may be able to confirm whether quasars reionize 
Hell at z sw 3. In addition, future observations of the 
HI and Hell Lya forest will soon be available with Sloan 
Digital Sky Survey III 25 and the Cosmic Origins Spectro- 
graph on the Hubble Space Telescope. These telescopes 
will significantly increase the number of sight-lines in the 
HI and Hell Lya forests. It is timely to make predictions 
for the effect of Hell reionization on these observations. 
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APPENDIX 

A. HEII ABSORBERS 

In this section, we examine the effect that dense systems have on filtering the radiation field and whether this effect 
is adequately captured in our simulations. We discuss two methods that we use to supplement the resolution of the 
simulations. 



Analytic Modeling 

The HI column density (Ahi) distribution is relatively well measured from the HI Lya forest, whereas the Hell 
column density (iVudi) distribution is significantly more uncertain. Given the well-measured ./Vhi distribution, we 
would like t o derive the A^Heii distrib ution and then use this to calculate the effect of filtering. We follow the method 
described in lHaardt fc Madaul ([1996). This approach uses the value of Ahi for an absorber and the incident values of 
Thi and Fueii to predict the value of iVncii- 

If the absorber is optically thin to 4Ry photons, it is trivial to infer the value of A^HcH, namely 



Audi = Ahi Jfthin, 



where 



??thii 



5-5 Fhc 
4 (1 - Yj 



He; 



1HI 

^HcII 



(Al) 
(A2) 
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and where Yn e is the mass fraction in helium and the factor 5.5 is the approximate ratio of Hell and HI recombination 
rates. HI column densities with 

JV H i,thin < 4 x a H l if 1 = 1.3 x 10 16 (—\ cm" 2 , (A3) 

yield optical depths of r < 1 for 4Ry photons, where n = Nn e u/Nm. 

What values of r\ do we expect? Taking th e T- values inferre d from the HI Lya and Hell Lya forests at z « 2.5 of 
Thi « lCT 12 s" 1 and r Hc ii w KT 14 s" 1 (e.g.. lShull et al.ll2004h . equation J22) yields ?7thin « 50. Or, ta king the z » 3 
value s of Thi ~ 10~ 12 s _1 and ThoII ~ 10~ 15 s _1 - at least in some regions - ^thin increases to 500 (e.g.. lReimers et all 
fl99l . 

An absorber is optically thick to 4 Ry photons while also being optically thin to 1 Ry photons for A^Hi,thin ^ -/Vhi ^ 
ctjjj = 1.6 x 10 17 cm~ 2 . Consequently, ThcII can be much lower within an absorber than the incident value, and, 
therefore, iVHeii can be la rger relative t o Nui than in the optically thin case. Thus, to calculate Nb_ q u we solve the 
quadratic formula given in lFardal et alJ (|1998f ) (their eqn. All), which determines A^jjeii assuming a slab of primordial 
gas. This equation requires as input -/Vhi, Thi, ThcII, a nd nn - the number dens ity of atomic and ionized hydrogen in 
the absorber. Th is method is an im provement over the lHaardt fc Madaul (|1996f ) multi-zoned model for incorporating 
Hell absorption (|Fardal et a l. 1998); 

To proceed with solving the lFardal et al.l (|1998h equation for -/VhcII, a model for n# (Ahi) is required. Simulations of 
the HI Lya forest have had remarkable success at matching the Nm distribution found in HI Lya forest observations 
fe.g.. lCeneraIlll994tlMiralda-Escudeet al.|[l99llHernauist et alj |1996; Ka tz et al.lll996l ), and analytic models have 
been constructed to understand these simulations. In particular, ISchavd ( 20011 ) argued that Lya forest absorbers 



have sizes that are comparable to the Jeans length, Lj(Sb,T). With the assumption that iVni = Hhi Lj and of 
photo-ionization equilibrium, it is trivial to map between tih and Ahi for JVhi ^ <?hi i namely (|Schavell2001h 



W \V3 / r \ 2/3 

-4 ,3 / iV HI \ r^O.17 I i HI 



"* « 6 x 10- c,^^^ ) T *[W^) ' (A4) 

where T4 — T/10 4 K. This formula agrees well with simulations (|Schavdl2001[ ). For our purposes it may reasonable to 
use equation IA4I even when ./Vhi ^ Ojjj , which is beyond its realm of applicability, because such high column density 
absorbers play an insignificant role in filtering the radiation near the Hell Lyman- limit. 

Figure [T751 plots the predictions of this model for r)(Nm) (bottom panel), the m.f.p. to be absorbed in systems 
with column density greater than ./Vhi (top panel), and the contribution of different systems to the m.f.p. (dA m f p (> 
Ani) _1 /rflog TVhi (middle panel). The function dA m f p (> iV"Hi) _1 /rf^VHi is proportional to the probabilit y a photon 
is abs orbed in a system of column density ./Vhi- These panels use the Al fit to the ./Vhi distribution in iFardal et al.l 
(1998). 26 All curves assume r H i = 10~ 12 s" 1 and T = 10 4 K. In the top panel, the solid curves take r He ii = 10" 15 s _1 , 
and the dotted curves take THeii = 10~ 14 s _1 . The vertical lines in both panels correspond to the approximate value 
of ./Vhi at which an absorber becomes self-shielded to Hell ionizing photons for ThoII = 10~ 14 s _1 (solid line) and 
THeii = 10~ 15 s _1 (dotted line). For the column densities at which an absorber becomes self-shielded, the value of r\ 
increases. However, this increase is compensated by the monotonically increasing nature of ?t.h(-/Vhi), and ultimately 
this effect wins out and J](N}n) decreases. This decrease is important, allowing 100 eV photons to have long m.f.p. s 
to be absorbed in dense systems. 

For r Hc ii = 10~ 14 s" 1 , absorbers with A H i > 10 16 cm~ 2 limit the m.f.p. of 4 Ry photons to «s 40 Mpc at the 
Hell Lyman-limit (see thickest, solid curve in middle panel in Fig. \W\\ . For ThoII = 10~ 15 s _1 , lower column density 
absorbers with 10 14 > JVhi ^ 10 15 cm -2 limit the m.f.p. of these photons to w 10 Mpc (dashed curves in middle and 
top panels). Higher energy photons have longer m.f.p. values. 

Filtering in our Simulations 

The radiative transfer grid employed in our study may not capture all of the absorbers with 7Vhi > 10 15 cm -2 
(Appendix D). These absorbers can play an important role in filtering the radiation field. Furthermore, even if these 
absorbers are resolved on the grid, the code will not properly handle the filtering when ThoU ^ Ai _1 . We describe 
below two methods that we use to better capture these effects. 

Filtering Method A 

The simpler of the two filtering methods uses the density structure on the radiative transfer grid and a modification of 
the radiative transfer code to achieve the appropriate filtering. Appendix |Pl shows that the column density distribution 
of absorbers measured from the 256 3 and 512 3 N-body simulation grids reproduce to within a factor of 2 the column 
density distribution seen in gas simulations. Therefore, our gridded N-body simulation will produce the correct filtering 
of the radiation field if it passes three tests: 

26 [Fardal °t al. ( 199Q) derived best fit functional forms for the distribution of absorbers as a function of ./Vhi an d z (see Table 1 in 
IFardal et al. 199SJ), fitting to a compilation of Lya forest observations. These fits account for the observed depression in d 2 N/dNnidz 
centered at Nm ss 10 16 cm~ 2 (Pctitiean et al. 1993). Note that the piecewise form for the curves in the middle panel in Figure [TrJl owes to 
the piecewise power-law that Fardal et alj C1998J) used to fit the observed Nm distribution. 
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Fig. 16. — Dotted curves in each panel are for Tho 
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and the solid are for ThoII 
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s x and all curves use 



Thi = 10 s . Top Panel: The m.f.p. for photons at select energies to be absorbed in a system with column density greater than -/Vhi- 
Middle Panel: The derivative with respect to Afjji of the inverse of the m.f.p. This quantity is maximized for the Nm that contribute most 
significantly to the absorption. Bottom Panel: The value of r\ = iVHell/^Hl in this model as a function of Nm. The vertical lines in both 
panels correspond to the TVhi at which an absorber becomes optically thick for the ThcII with the same line style. 



The density structure of the absorbers is captured to the level that is required. 

The radiative transfer method accounts for the filtering of the radiation field as radiation traverses overdense 
cells. 



• The clustering of the absorbers is taken into account. 

When the absorber is optically thin, r\ is independent of the absorber density. As seen in the middle panel in Figure [TBI 
a significant fraction of the filtering owes to optically thin absorbers when ThcH = 10~ 15 s _1 . Even when an absorber 
becomes optic ally thick, the dependence of 77 on density is weak for the most relevant column densities if one uses the 
ISchavd (|2001h model for tih(^Vhi)- (Note the flatness of r\ in the bottom panel in Fig. [16] for the JVm that contribute 
the most to the m.f.p.) Therefore, since the iVm column density distribution is captured fairly well at our grid scale, 
there is reason to believe that the density is adequately described for the purpose of filtering the radiation field. 

Passing the second of the three tests requires an addition to the radiative transfer algorithm beyond what was 
outlined in the previous section. For r^ eII <C At, a cell in the code outlined in [J2] would have XHcii < ^Hcii.cq once 
some fraction of the rays that will travel through that cell in a timestep have done so. Therefore, the filtering will 
be underestimated for the rays that enter after xucii < ^Hcii,cq- To remedy this, once a cell with A5 > Ap has 

£Hcii < ^Hcii.cq, we perform the radiative transfer as if xhcII = ^Heii.eqj where XHeii,eq = Qa n e /Tn e ii and Fneii 
is tabulated from all the previous rays that have entered the cell within the timestep up to the current ray. Note 
that £Heii always approaches £Heii,eq from below because recombinations are performed prior to ionizations within a 
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timestep. We set Ap = 4, which for Thi = 10~ 12 s -1 corresponds to 



W „, = 2 x 10-W ('¥)' (VgbH 1 ^ h (A5) 



'cell \ 

IMpc/ 

Passing the final of the three tests - capturing the clustering of the absorbers - is the most difficult and is not 
achieved with this method. Even though the column density distribution is maintained, the bias of the absorbers is 
changed by the coarse gridding (Appendix D). ft is conceivable that capturing the bias is not of crucial importance 
because the m.f.p. for ionizing photons is typically larger than the correlation length of the absorbing systems. This 
method does not filter the background radiation (> 500 eV). Figure [TBI suggests that such filtering is not important 
because the m.f.p. for these photons to be absorbed in a dense system is gigaparsecs. 

Filtering Method B 

The secon d method fo r filter ing the radiation is based on the model proposed in lHaardt &: Madaul (|1996f ) and 
improved in Farda l et al.l (|f998f ) to infer -/VhcII from ATjji, which is described in the beginning of this section. This 
inference also requires as input ThoII- We take it to be the maximum of the current tabulated value of ThoII m the 
cell and the value of THeii from the previous timestep. 

In addition, this method assumes that Lya absorption arises from gas at the outskirts of dark matter halos. We 
assume that the cross section for a photon to intersect an absorber of column density A^eii associated with a halo of 
mass m has the functional form 

a(m, N-aai) = o m , P(N HoU ) (™/m») 2/3 , (A6) 

where -P(AhcIi) is taken to be the probability distribution of iVneii that corresponds to lines with log 10 (iVHi) > 
14.5 and is zero for smaller column densities, and o~ mt sets the normalization. In this method, rays that have 
photons of energy E 1 and that travel through a cell with a halo of mass to are attenuated by the factor 
/ dN HoU a(m, N HcU ) exp[-r(£ 7 , iV He ii, AWi)]/^, where r = ct H i(£ 7 ) N m + a" He ii(-E 7 ) N HcU . 

The to 2 / 3 scaling in equation (|A6j) is chosen such that the cross section is proportional to the square of the 
halo virial radius. The values of a m , and to* are assigned so that the mass integral over cr(m, iVneii) n h{m), 
where nh(m) is the halo mass function, yields dN / dldN^ c u - the number of abso rbers per proper le ngth I per 
Audi- To derive dN/dldN Rcn , we use the Al fit to dN/dldN m (z) provided in iFardal et all (|1998t) and that 
dN J dldNncii(z) = ?7(-/Vhi) _1 dN / dldNm{z) . While we use the observations to construct -P(AhoIi), physically it has 
to do with the density profile around halos. The motivation for this filtering method is tha t the high column density 
absorbers are shown to be associated with halos (|Cen fc Simcodll997l : lMiralda-Escudelfl998f ) . This method filters both 
the rays and the background radiation. 

B. TESTS OF CODE 
Number of Photon Bins 

Since the computation time scales linearly with the number of frequency bins n v and since memory requirements 
also increase linearly with n v , it is important to find the minimum value of n v that achieves acceptable convergence. 
To test how our code converges with the number of frequency bins, we place a quasar in a homogeneous f GM initialized 
with XHoii = f and solve for XHeii and T while varying n v . Figure [T71 displays the temperature (top panel) and ^Hein 
(bottom panel) for n v — 3, 5, and 50. For simplicity, we have set the speed of light to infinity. The left-most curves 
are for a source with auv = 1.5 and N = 10 54 ionizing photons s _1 after t — 50 Myr and the right-most curves are 
the same after t = 250 Myr, where t — corresponds to z = 4. The curves at fixed t have only minor differences 
on scales that have xhcIH ^ 0.01. However, for xhcIH < 0.01, 3 frequency bins provides fairly poor convergence to 
the true solution (which essentially is represented by the n v = 50 case), whereas 5 frequency bins shows much better 
convergence. The cosmological simulations presented in this paper use n v — 5 unless otherwise specified. 

Even though for a single source the solution is well converged to fairly large radii for n v — 5, in the simulations the 
small inaccuracies at large radii add up and can result in a significant error. Figure [18] compares at Sfjoiii = 0.85 of 
the fiducial 5 frequency bin simulation - simulation Dl - to one with 10 frequency bins - D3 - with the bins centered 
at (57, 66, 76, 90, 107, 129, 162, 209, 283, 409) eV such that each bin carries an equal energy. We have checked that the 
latter simulation is well converged to the full solution by running simulations between 5 and 10 bins and with different 
spacings in energy. The ionization field is similar between simulations Dl and D3, with the voids being slightly more 
ionized in Dl. The value of T is 300 K greater in Dl, and it can be seen that some of the hotter regions are a bit 
hotter in Dl than in D3. However, these differences are small compared to the differences that arise owing to, for 
example, the uncertainty in auv- We have also compared the temperature PDF between Dl and D3, and the PDFs 
are almost identical. The temperature PDF still extends to 35, 000 K in D3 as in Dl at the end of Hell reionization. 

Code Comparison 

Analytic solutions to multi-frequency radiative transfer problems generall y do not exist. To test our code, we 
compare it with the 1-D radiative transfer code presented in lLidz et al.l (|2007l ). To simplify the comparison, we use a 
1-D version of our code that is identical t o our 3-D code exc ept that it sends a single ray during a timestep to calculate 
the spherically symmetric solution. The iLidz et all (|2007[ ) code was written independently by Adam Lidz who was 
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Fig. 17. — Test of convergence with the number of frequency channels for the case of a single QSO in a homogeneous IGM with xjjell = 1. 
The top panel shows the temperature, and the bottom shows the Helll fraction. The left-most curves are for t = 50 Myr after a source 
with qtjv = 1-5 and N = 10 54 ionizing photons s — 1 turns on and the right-most curves are t = 250 Myr after. 
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Fig. 18. — Comparison of the fiducial 5 frequency bin simulation 
panels are taken from the z = 3.3 snapshot, for which xhoIII = 0.85. 



simulation Dl - to a simulation with 10 frequency bins — D3. The 



not involved in the development of our new code. The methodology used in our code is significantly different from 
that used in the lLidz et al.l (|2007f ) code. The latter code assumes a power-law spectrum, and uses this power-law and 
the intervening column densities of HI, Hel, and Hell to infer the incident spectrum at some radius from the source. 
Furthermore, it uses much smaller timesteps than our code and iterates to achieve convergence. 

We analyze the solutions provided by these codes for the simple case of a quasar that produces 10 54 s _1 Hell ionizing 
photons with a spectral index of 1.5 in a homogeneous IGM with XHeii = 1- We solve for the ionization front of this 
quasar assuming a homogeneous, static (non-expanding) universe, with density equal to the mean density at z = 4. 
Initially, the hydrogen is assumed to be ionized and the helium to be singly ionized. For this calculation, the speed 
of light is taken to be infinite, the code used in this paper is set to have 1 Myr timesteps, and the same ray-casting 
parameters as used in the cosmological runs presented in this paper. 

Figure \T§\ compares the evolution of t he temperature an d SHeiii front as a function of time between a 1-D version of 
our code (dotted curves) and the code in lLidz et al.l (|2007l ) (solid curves). The left-most curves are the temperature, T, 
and ShcIII values as a function of radius after 2 Myr (2 timesteps for our code), the middle curves are those after 8 Myr 
(8 timesteps), and the right-most curves are those after 20 Myr (20 timesteps). At r < 5 Mpc, the agreement between 
the T{r) produced from the two codes is not perfect (top panel). The slight disagreement owes to a simplification our 
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Fig. 19. — Comparison of the temperature (top panel) and Hell ionized fraction (bottom panel) between our code (dotted curves) and 
the code in[Lidz ct al. (2007) (solid curves). The QSO is shining in a homogeneous medium at z = 4 with N = 10 54 Hell ionizing photons 
s — . The left-most curves are T{r) and ShcIIiM a fter 2 Myr, the middle curves are these functions after 8 Myr, and the right- most curves 
are after 20 Myr. 

algorithm makes: it groups all of the photons released in a timestep into a single photon bundle. In reality, harder 
photons will travel ahead of softer photons because softer photons are a bsorbed first. We have verified that if we 
decrease our timestep the solution from our code converges to that of the iLidz et ail (|2007l ) code. However, our 3-D 
code sends multiple rays for each cell at the radii where the discrepancy exists, capturing better the region where 
the 1-D code has not converged. We have verified that this minor discrepancy essentially disappears in the 3-D code. 
Furthermore, even though timesteps are larger in the cosmological runs than in this test, the convergence of this code 
is most dependent on the number of timesteps rather than the timestep duration or source luminosity. 

The bottom panel in Figure [T9l compares the Helll ionization fronts in the two codes. We plot the y-axis in log 
because the solutions for the fronts are essentially indistinguishable on linear scales. At large radii, there is a minor 
systematic difference in the ionization level. This minor difference is much smaller than other uncertainties inherent 
in Hell reionization calculations, in particular the value of afjv- This systematic does not depend on t he number 
of frequency bins (for n v > 5) or decreasing the timestep in our code. Also, note that the output of the ILidz et all 
(2007) curve is at 1.9 Myr rather than 2 Myr for its left-most curve, which accounts for some of the discrepancy. 
The maxim um differences between the 8 Myr curves is 3%. The small differences at large r could arise because the 
ILidz et all (|2007f) code follow the full three level system, whereas our code ignores Hel. 

Temperature- Density Relation 

Here, we investigate how well our scheme for computing the temperature evolution in N-body simulations works. 
The most difficult aspect of such a scheme is to evaluate the convective derivatives that account for the flow of matter 
through cells. Our scheme for capturing this is discussed in Section 12.21 To test our code, we compare the resulting 
T — At, relation in a simulation without Hell reionization to that of the iHui fc Gnedinl (|1997t ) analytic model. This 
analytic model has been shown to reasonably match the relation found in semi-analytic calculations that use the 
Zel'dovich app roximation and the e volution of the temperature density relation found in hydrodynamic simulations 
(although, the lHui fc Gnedinl (|1997[ ) comparison is rather limited for the latter case). 

Figure 1^01 compares the evolution of the T — A;, equati on of state 7 — 1 in o ur code using the 256 3 grid in the 186 
Mpc box to the 7 — 1 predicted by the analytic formula in lHui fc Gnedinl IJ1997I ) for models where z rcl = 6.2. Our 7 — 1 
are calculated by fitting the T — A^ relation to a power-law, weighting the fit by the number of grid cells at each A&. 
We have checked tha t our results do not change if we instead fit only to A& values near unity, which is closer to what 
IHui fc Gnedml (|1997ft calculate. 

The agreement is acceptable between the two predictions for 7 — 1. We would not exp ect exact agreement. Our code 
calculates the temperature evolution on a non-linear scale (the cell scale) whereas th elHui fc Gnedinl (1997) formula 
uses linear theory. Therefore, our value for 7 — 1 should evolve more quickly than the IHui fe Gnedinl (|l997l ) estimate 
since A& on average grows faster on non-linear scales. 

Our scheme misses the effects of shocks, which add dispersion to the T — A& relation. The s.d. in the T — A& relation 
fro m shocks is < 500 K at Ab = 0.5 and < 1000 K at Af, = 1, and it increases steadily as Af, increases (see Figure 
1 in lHui fc Gnedinl 1 19971 ). The dispersion in temperature at low densities in simulations that include shocks (and a 
uniform ionizing background) is much smaller than the dispersion that results from Hell reionization. 

C. IONIZATION CORRECTION 
The fraction of helium in Hell evolves according to the equation (assuming Hel is in ionization equilibrium) : 

d^Hell 



dt 



FhcII #HeII + a A n e XHcIII- 
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Fig. 20. — Evolution of 7 — 1 in our code compared to that predicted by IHui fe Gnedinl (|1997T ) (HG97) for T rci 
T rei = 10, 000 K. 
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Fig. 21. — Error in our code's calculation of xn c u(At) in units of the number of recombinations plotted as a function of THell- This 
calculation assumes At = 20 Myr. The blue line represents the case £Hell(0) = 0.01 and S b = 0, the green line is for £Hcll(0) = 1 and 
(5;, = 0, and the red line that for £hcIi(0) = 1 and <5f, = 10. 

If we take n e to be constant (a relatively good approximation during z ~ 3 Hell reionization because hydrogen is 
ionized), equation (IClj) yields 

EHeIl(*) = ZHcII.cq + (»HeIl(0) - ^Hcll.oq) exp [-t/t cq \ , (C2) 

where t oq = (Fh c ii + l/t r0 c) _1 , ^Heii(O) is the initial Hell fraction, and 



^Hell.i 



a A (T) n e 



Hell 



a A (T) n e 



(C3) 



Compare the exact solution for XHeii(t) given by equation (|C2|) to that without recombinations (and where an 
estimate for the number of recombinations qa n e XHcii(O) At is added to XHeii at the beginning of the timestep): 



^Heii(i) = [jCHeii(O) + a A n e XHoiii(O)] exp [-r Ho ii t] . 



(C4) 



The value XHoii(i) is what is computed by our Hell reionization code over a timestep At. The computation of xn e u(t) 
allows us to determine the effect of all sources on a single cell independently - ThoII in equation (|C4p is the sum of 
the ThoII from all the incident rays - aside from the issue of the order rays from different sources reach a cell which 
affects the heating rate in the cell and how the cell shadows other cells. 

Figure [21] shows that the error generated by using equation (|C4[) rather than the true solution, equation (|C3|) , is 
minor. This figure plots the ratio of Accneii to an estimate for the number of recombinations in a timestep «a n e At 
as a function of Tneii for At = 20 Myr. The blue line is for XHcii(O) = 0.01 and Sb = 0, the green line for XHcii(O) = 1 
and Si, = 0, and the red line for XHcii(O) = 1 an d Si, = 10. The error in Axneii m all cases is small, less than 0.3 in 
units of the number of recombinations (and AxhcII ^ 0.005). For other values of At, we find that AxHeii/(o ; A n e At) 
peaks at roughly the same amplitude of 0.3, but with the peak at ThoU ~ 1.8 Ai~ . 

This method also reproduces the correct heating rate with minimal error. It is easy to check this in the case of 
ionization equilibrium, the regime where one might suspect that this code is unreliable. In ionization equilibrium, the 
amount of heating is 

AQ - r Hc ii (AE) XHdi.cqAt = a A n e (AE) At, (C5) 
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Fig. 22. — Comparison of the dark matter and gas volume-averaged PDFs for z = 3.5. The solid curves are the dark matter PDFs and 
the dashed curves are the gas PDFs measured in the specified cell sizes. The dot-dashed curves are the MHR gas PDF. 

where (AE) is the average excess energy above 4 Ry of a photon that enters a cell and i oniz ed a Hell ion, and we have 
used equation IC3I for XHeii.eq hi the limit ThoII 3> ctA^E- Comparing AQ in equation IC5I to the value our algorithm 
provides 

AQ = r H cii (AS) Adduce = a A n e {AE) At, (C6) 

where we have used that AxHcii,rcc ~ Q^A.Hoii n e At (since £Hcii ~ 1) and that all the Helll ions that recombine are 
reionized (i.e., r H oii > lArcc)- 

Since an accurate estimate for xhcII is also necessary for Hell Lya forest calculations, we correct the mis- 
take our approximate method incurs by substituting the true solution, i.e. equation JUTJ, for cells that have 
THeii,Lya(^Heii, z, <5&) > 10 (which correspond to larger values of THcii.Lya^Hcii, z, St)), using the formula 



THcU,Lya{5:Hcn,Z,Sb) — 3.4 



£HcII 

IF 3 



1 



3/2 



A,, 



(C7) 



This correction might be worrisome because the algorithm no longer consistently calculates the balance between 
recombinations and ionizations. However, a minor correction of AxhcII = I^Hoii(Ai) — iHeii(At)| in a cell is acceptable 
because it is smaller than other uncertainties, such as uncertainties in the recombination coefficient (Section 12.5ft . 

While we have used Case A recombination coefficients in the tests above, the equivalent discussion applies if we replace 
Case A with Case B, as is normally done by our code (except Axneii will be smaller for Case B). We always use the Case 
A coefficient to correct our code's calculation because ground state recombination photons are not absorbed locally 
in regions with significant transmission. Rather, these photons are absorbed in dense (probably neutral) systems. 
As a consequence, our calculation of THeii.Lya misses the contribution to T from ground state recombinations and, 
therefore, underestimates the amount of transmission. We can put an upper limit on the contribution to ThoU from 
such recombinations. If we assume that the volume-averaged photo-ionization rate is just large enough to balance 
the number of recombinations aside from those to the ground state and that 55 eV recombination photons have the 
same m.f.p. as the typical ionizing photons from quasars, then r reC: Hdi = ("A — «B)/aA IhcII ~ 0.3rHcii- In °ur 
simulations, the photo-ionization rate is larger than what is required to balance the number of recombinations (except 
in select regions) and the m.f.p. for Hell-ionizing recombination photons is shorter than for other photons, such that 
r rC c,Hcii < 0.3rHeii- 

D. DARK MATTER VERSUS GAS DISTRIBUTION AND THE N m DISTRIBUTION 

The radiative transfer computations in this study are performed as post processing on top of gridded N-body fields. 
This assumes that the gas fluctuations are pressure smoothed at roughly the simulation grid scale and trace the dark 
matter on larger scales. Cold dark matter clumps on extremely small scales, and, therefore, it is important to properly 
choose the scale at which to grid the N-body field to approximate the gas pressure smoothing. Here, we justify our 
grid sizes. 

We use the Q6 simulation of ISpringel fc Hernquistl (|2003f h This simulation consists of a 14.3 Mpc box with 486 3 
dar k matter particles and i nitially 486 3 SPH particles. The Q 6 simulation includes gas cooling, s tar formation, 
the lHaardt &; Madaul (|1996( l m odel for the ionizing backg round IjKatz et al.l Il996ak iDave et all [l99S) , as well as a 
prescription for galactic winds (jSpringel &; HernquistJ 120021 ). It resolves the Jeans scale at mean density with rs 10 4 
SPH particles, and it has a force resolution of 1.2 proper kpc. The following analysis uses the z = 3.5 snapshot from 
this simulation. 

We first calculate the dark matter and gas overdensity PDFs to quantify how well the gridded dark matter traces 
the gas. This is shown in Figure |2"21 for 0.45 and 0.89 comoving Mpc grid cells. These cell sizes are comparable to 
the respective cell sizes for the 512 3 and 256 3 grids in the 186 Mpc simulation box. The solid curves are the dark 
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Fig. 23. — Number of absorbing columns per z per log^A^i- The blue dotted curves represent the case in which this quantity is 
measured in 0.89 Mpc coarse cells. The thicker curve uses the average dark matter density in the coarse cell to measure Nm, and the thin 
curve uses the gas density along a skewer the length of the coarse cell of the pressure-smoothed gas. The green dashed and red solid curves 
are the same as the blue dotted but measured in 0.45 and 0.22 Mpc coarse cells, respectively. 



matte r PDFs on the grid and the dotted are the same for the gas. The dot-dashed curve is the iMiralda-Escude et al.l 
(2000) (henceforth MHR) gas PDF. This PDF represents that of the pressure smoothed gas (grid artifacts should be 
unimportant) . The extent to which the gridded dark matter PDF agrees with the MHR gas PDF reflects how well the 
dumpiness of the gas is captured in the simulations. 27 The low density gas PDF is important for accurately calculating 
from the simulations the absorption in the Hell Lya forest, the number of recombinations per Hell-ionizing photon, 
as well as the filtering of the radiation field. Figure [22l demonstrates that the gas dumpiness is described reasonably 
well with the chosen grid scales. 

The high-density density regions, in particular densities responsible for absorbers with Ar-i > 10 15 cm -2 , are im- 
portant for filtering the Hell radiation field (as shown in Appendix A). Figure [2"o1 compares the Ahi column density 
distribution of the coarse cell-smoothed dark matter (thin curves) to the column density distribution of the pressure- 
smoothed gas (thick curves) , with both column densities measured along skewers that have the same length as a coarse 
cell. Note that Am = a;Hi, C q 'cell nii(z) A&, where we use Thi = 10~ 12 s _1 . The blue dotted curves represent the 
case in which Ar-i is measured with 0.9 Mpc coarse cells. The green dashed and red solid curves are the same but 
measured in 0.45 and 0.22 Mpc coarse cells, respectively. The dark matter field, which for coarser grids underpredicts 
the number of Ar-i > 10 15 cm~ 2 absorbers, overproduces the abundance of absorbers for 0.22 Mpc cells. Interestingly, 
the Ahi distribution of the 0.45 Mpc gridded field (which corresponds roughly to the 512 3 grid simulation in the 186 
Mpc box) agrees well with the gas distribution. 

E. QSO MODELLING 

QSO Spectral Index 

Both methods for populating the box with QSOs discussed in Section [3] return the specific luminosity at the HI 
Lyman-limit. We extrapolate from this spe cific luminos i ty to higher energies using a power-law with index auv- 
The parameter any has been measured by iTelfer et al.l (|2002D to be ~ 1 .6. but with an approximately Gaussian 
distribution with a large standard deviation (s.d.) of 0.86 (|Telfer et al.ll2002n. The large s. d. in auv amon g QSOs may 
help explain transmission fluctuations in the Hell Lya forest at z w 2.5 (Sh ull et al.ll20 04). More recently. IScott et al.l 
(2004) measured a much smaller mean auv from a sample of both FUSE and H S T quasars. They d erived an average 
value of auv ~ 0.6 and with a variance compara ble to that o f Telfer et ail (|2002t ). IScott et ail (2004) also find a strong 
correlation betw een auv and QSO luminosity. IScott et ail (|2004h argued that the discrepancy between their mean 
auv a nd that of ITelfer et al.l (|2002[) likely owed t o their sample containing lower luminosity quasars. 



The ITelfer et alJ(|2002D and IScott et all (|2004[ ) values of auv are derived from fitting 1 Ry < E 1 < 4 Ry. For 



Hell reionizatio n, the spectr u m be tween 4 Ry and 1 ke V is relevant . If we naively extrapolate from 1 Ry to the soft 
X-ray using the ITelfer et al.l (|2002l ). or particularly the IScott et al.l (|2004l ). auv distribution this would result in an 
over-estimate of the soft X- ray luminos i ty fun ction for quasars because of the large variance these studies measure 
in auv (and in the case of IScott et al.1 (|2004h . the low value of auv)- The average spectral index ao tX needed to 
j oin the luminosity at 2500 A and 2 keV has been meas ured from a sample of optically selected quasars to be ~ 1.6 
(Stcff en et al.l l2"006L 28 Furthermore, Steffcn ct al. ( 200 a) m easu red a dispersion in this index 0.08 — 0.14, significantly 
smaller than the dispersion in auv that Tel fer et al.l (2002) and IScott et al.l (|2004f ) infer. 



27 Note that the MHR PDF is generated using Eulerian simulations that employ a different cosmology and thermal history from simulation 
Q6. While the gas PDF is fairly robust to the cosmology, the corresponding PDF of the pressure smoothed gas in simulation G4 will be 
slightly different than the MHR PDF. We choose to use the MHR PDF because we find that particle noise in SPH simulations makes it 
difficult to capture the low density PDF. 



' Stcffen ct al. (2006) found that this power-law was correlated with luminosity, ranging from 
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Escape Fraction 



Not all of the ionizing radiation escapes from a quasar. The fraction of Hell-ionizing radiation that is obscured is 
denoted by f cov . Measureme nts in the X-ray s uggest f cov rj 0.5, with evidence that / cov increases with luminosity 
(|Gilli et al.ll2007t ). We use the iGilli et all (|2007f ) fitting function for / cov in calculations that require it. 29 

The fraction of ionizing radiation that escapes from QSOs is almost certainly directionally dependent - it will be 
easiest for ionizing photons to escape on trajectories that do not intersect infalling material or galactic disks. In select 
simulations, we take a beam for the ionizing radiation from the QSO with solid angle 51b = 27r (1 — f cov ) along both 
axes. We assume the symmetry axis of the beam is randomly oriented. For the other simulations, we uniformly 
suppress the intensity along all sight-lines from the QSO by the factor / CO v This approximates the scenario in which 
the obscuration owes to infalling material and (when averaged over the QSO lifetime) is isotropic. The value of the 
covering fraction is only required for Method I if the quasars are beamed. It affects the quasars in Method II by 
suppressing their ionizing luminosity. 
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